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ABSTRACT 


A  numerical  study  of  a  generic  axisymmetric  ramjet  operating  at  conditions 
corresponding  to  flight  Mach  3.0  and  a  standard  altitude  of  10  km  is  presented.  The  study 
includes  both  modeling  of  steady-state  flowfields  in  the  ramjet  as  well  as  transient 
throttling  maneuvers  in  which  the  throttle  is  decreased  or  increased  from  maximum  or 
minimum  throttle  positions.  The  results  presented  here  focus  on  entropy  generation  and 
performance  characteristics.  Combustion-generated  exothermic  heat  release  is  modeled 
using  simple  volumetric  energy  addition  to  the  flow  within  a  defined  heat  release  zone. 
The  study  utilizes  two  levels  of  wall  boundary  modeling,  corresponding  respectively  to 
inviscid  and  viscous  walls  in  the  ramjet.  The  second  level  of  modeling  (with  viscous 
walls)  presents  many  challenges  due  to  the  inherent  tendency  of  the  no-slip  boundary 
condition  to  cause  reverse  flow  to  develop  in  the  ramjet,  particularly  along  the  wetted 
surfaces  of  the  inlet  where  the  adverse  pressure  gradient  associated  with  the  deceleration 
and  heat  release  in  the  ramjet  has  the  largest  initial  impact.  This  separated  flow  results  in 
eventual  unstart  of  the  ramjet  due  to  large-scale  propagation  of  the  separation  upstream; 
there  is  also  inherent  unsteadiness  due  to  boundary  layer  effects. 

To  address  the  challenges  presented  by  the  no-slip  boundary  condition,  a  bleed 
boundary  condition  specified  at  the  inlet  throat  is  incorporated.  This  bleed  extracts 
approximately  10%  of  the  mass  flow.  As  an  alternative  to  bleeding  mass  from  the  flow 
path  of  the  ramjet,  a  generic  (alternative)  model  of  a  ramjet  dump  combustor  is  also 
studied.  This  configuration  has  a  geometry  in  which  a  constant  area  heat  addition  zone  is 
located  downstream  of  a  large  step  at  the  exit  of  the  ramjet  inlet.  This  configuration  is 
analyzed  and  compared  to  the  original  configuration  without  the  dump  combustor.  It  is 
found  that  both  the  bleed  boundary  condition  and  the  dump  combustor  are  extremely 
effective  at  preventing  the  normal  shock  from  propagating  upstream. 
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NOMENCLATURE 


Symbol  Description 


SHt  Entropy  Generation  Per  Mass  Due  to  Ideal  Heat  Transfer 

5Q  Amount  of  Uniform  Heat  Rate  Addition 

Ttot  Total  Temperature 

SnHT  Entropy  Generation  Per  Mass  Due  to  Non-Ideal  Heat  Transfer 

T  Ambient  Temperature 

Ax  Local  Area  at  Station  "x" 

Ax  Chocked  Flow  Throat  Area  at  Station  "x" 

y  Local  Gas  Constant 

Mx  Local  Mach  Number  at  Station  "x" 

px  Local  Pressure  at  Station  "x" 

p0x  Local  Total  Pressure  at  Station  "x" 

T0x  Local  Total  Temperature  at  Station  "x" 

m  Mass  Flow  Rate  Entering  Ramjet 

CP  Constant  Pressure  Specific  Heat 

Q  Heat  Addition  Corresponding  to  Maximum/Minimum  Throttle  Settings 

HT  Heat  Transfer 


1.  INTRODUCTION 


Accurate  transient  analysis  and  early  incorporation  of  the  effects  of  transients  in 
the  design  and  optimization  of  high  performance  military  aircraft  are  increasing  in 
importance  due  to  the  demanding  nature  of  proposed  aircraft  missions  and  the 
development  of  novel/non-traditional  weapons  systems  [1].  The  effects  of  both  proposed 
missions  and  new  weapons  systems  are  tied  directly  to  the  aircraft  propulsion  sub-system 
in  terms  of  the  engine  providing  required  performance  margins  and  also  to  the  strong 
fluid  and  thermodynamic  interaction  and  integration  of  the  propulsion  sub-system  with 
other  vehicle  sub-systems  [2].  The  growing  complexity  of  proposed  aircraft  and  the 
pressing  thermal  challenges  presented  by  defined  mission  requirements  demands  aircraft 
optimization.  Optimization  for  both  preliminary  and  final  design  is  increasingly  being 
done  via  multi-disciplinary  analysis  and  optimization  (MDA/MDO)  [3].  MDA/MDO 
allows  the  entire  aircraft  (often  as  it  is  flown  through  a  mission)  to  be  analyzed  and 
optimized  as  a  whole,  in  contrast  to  the  traditional  approach  of  optimizing  each 
component  and  sub-system  on  the  aircraft  and  then  integrating  the  separate  subsystems  of 
the  aircraft  together  (hence  at  least  nominally  creating  an  overall  suboptimal  aircraft 
system).  In  order  to  properly  allocate  time  and  money  in  areas  that  need  improvement 
when  optimizing  an  entire  aircraft  using  MDA/MDO,  a  common  loss  metric  would  be  of 
great  value  when  evaluating  the  vehicle-driven  inefficiencies  of  all  subsystems  across  an 
aircraft.  This  common  loss  metric  is  provided  by  analyzing  exergy  destruction  or  entropy 
generation  [4]  [5];  exergy  destruction,  when  correctly  done,  allows  all  subsystems  to  be 
analyzed  using  the  same  fundamental  metric  161 .  This  also  allows  for  more  effective 
allocation  of  time  and  resources  within  an  aircraft  design  and  optimization  process,  in 
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terms  of  identifying  where  losses  occur  and  the  actual  magnitude  of  their  effect  on 
vehicle  performance.  The  central  purpose  of  the  current  research  is  to  assist  in 
understanding  and  formalizing  the  analysis  of  entropy  generation  within  a  transient  event 
for  a  simple  generic  low- speed  engine  flow-field. 

The  overall  objective  of  this  investigation  is  to  examine  entropy  generation  in  a 
simplified  and  highly  generic  ramjet  model  with  focus  on  both  steady-state  results  and 
throttling/transient  responses.  This  work  will  serve  as  part  of  a  platform  for  future  work 
using  entropy  as  an  optimization  metric  utilized  in  MDA/MDO  approaches.  The 
objective  was  accomplished  by  utilizing  computational  fluid  dynamics  (CFD)  in  order  to 
analyze  the  fluid  dynamics,  entropy  generation,  and  performance  characteristics  for  three 
distinct  cases  of  ramjet  throttling.  The  first  case  examined  is  an  inviscid  wall  analysis  of  a 
conventional  axisymmetric  ramjet  design  with  a  faired  continuous  combustor,  i.e.  no¬ 
dump  or  step,  between  inlet  and  combustor.  The  second  case  is  also  an  inviscid  wall 
analysis  of  an  axisymmetric  ramjet,  however,  with  a  large  dump  combustor  [7]  [8].  The 
third  case  provides  results  of  a  viscous  wall  analysis  of  the  conventional  axisymmetric 
ramjet  configuration  used  in  the  first  case.  A  fourth  case  would  have  included  the  viscous 
wall  analysis  of  the  axisymmetric  ramjet  with  a  dump  combustor  (same  configuration 
modeled  in  the  second  case);  however,  the  results  of  the  CFD  analysis  proved  to  be  too 
computationally  expensive  and  the  data  would  not  have  been  available  within  a  timely 
fashion. 

Entropy  generation  per  mass  (simply  referred  to  as  entropy  generation  in  this 
work)  is  a  main  focus  in  this  investigation.  It  is  used  here  as  a  means  of  tracking  losses  in 
the  ramjet  flow-field.  Entropy  generation  is  allocated  in  terms  of  both  spatial  location  and 
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loss  mechanism.  The  latter  (physical  mechanism)  allocation  includes  entropy  generation 
due  to  friction,  convective  heat  transfer,  non-ideal  heat  transfer,  ideal  heat  transfer,  and 
shocks.  This  allocation  of  entropy  generation  into  its  components  allows  for  losses  to  be 
allotted  to  their  respective  physical  mechanism,  which  should  be  useful  in  MDA/MDO 
efforts.  Work  which  describes  in  some  detail  the  fundamentals  and  applications  of 
entropy  generation  via  distinct  physical  mechanisms  and  processes  has  been  previously 
done  [4].  The  contribution  to  entropy  generation  due  to  friction  is  the  result  of  shear  stress 
(from  velocity  gradients)  across  adjacent  stream-tubes.  Similarly,  entropy  due  to 
convective  heat  transfer  arises  from  the  thermal  gradients  across  adjacent  stream-tubes. 
The  entropy  associated  with  the  heat  addition  was  divided  into  two  categories,  ideal  and 
non-ideal  heat  transfer.  Entropy  due  to  ideal  heat  transfer  is  calculated  using  Equation 
(1). 


If  heat  were  added  to  a  stagnate  flow  (velocity  being  zero),  the  local  temperature  would 
equal  the  total  temperature  and  thus  all  of  the  entropy  generated  due  to  heat  addition  is 
considered  ideal.  However,  because  the  velocity  is  not  zero  across  the  heat  addition 
region,  a  component  of  entropy  arises  do  to  the  non-ideal  nature  of  heat  addition 
calculated  using  Equation  (2). 

S-ht=8Q  (t~)  <2> 

Observation  of  Equations  (1)  and  (2)  reveals  the  trade-off  between  ideal  and  non-ideal 
entropy  generation;  as  velocity  increases  (total  temperature  decreases  and  ambient 
temperature  increases)  the  entropy  associated  with  non-ideal  heat  transfer  rises  and 
entropy  due  to  ideal  heat  transfer  declines.  Note  that  if  the  total  temperature  were  equal  to 
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the  ambient  temperature  (stagnate  flow),  Equation  (2)  would  equal  zero  and  entropy 
generated  due  to  heat  addition  would  entirely  consist  of  entropy  due  to  ideal  heat  transfer. 
Entropy  due  to  shocks,  being  the  only  mechanism  remaining  contributing  to  entropy 
generation,  accounts  for  the  difference  between  the  total  entropy  and  the  sum  the  entropy 
due  to  friction,  convective  heat  transfer,  ideal  heat  transfer  and  non-ideal  heat  transfer. 

A  ramjet  was  chosen  for  this  analysis  due  to  the  inherent  simplicity  of  the  design 
(no  turbo  machinery)  and  also  because  of  their  relatively  large  flight  Mach  number 
ranges  [9].  Ordinarily,  ramjets  use  the  ram  effect  to  compress  incoming  supersonic  air  in 
preparation  for  combustion  by  use  of  a  converging  nozzle  [10].  The  flow  is  then  passed 
through  a  normal  shock,  or  a  shock  system,  which  is  located  between  the  inlet  throat  and 
the  main  combustor,  such  that  subsonic  flow  is  established  at  the  entrance  of  the 
combustor.  Because  of  high  temperatures  and  pressures  associated  with  the  ram  effect 
(and  the  presence  of  a  normal  shock),  mechanical  compression  via  a  compressor  or  a  fan 
(as  found  in  turbojets  and  turbofan  engines)  is  not  needed.  Thus,  a  downstream  turbine  is 
not  necessary  in  order  to  extract  energy  from  the  flow.  However,  for  this  reason,  ramjets 
operate  most  effectively  in  supersonic  flows  up  to  about  a  flight  Mach  of  5  or  6,  when 
thermal  issues  overwhelm  the  cycle.  This  reliance  on  the  ram  effect  prevents  ramjets 
from  being  used  as  a  propulsive  device  to  achieve  flight  without  high  speed  flow  entering 
the  inlet,  i.e.  ramjets  cannot  be  used  for  takeoff.  Several  methods  have  been  devised  to 
alleviate  this  problem.  Solid  rocket  motors  have  been  developed  that  utilize  the  ramjet 
configuration  as  storage  for  the  propellant  material  [9].  As  the  propellant  is  consumed  and 
flight  velocity  is  increased,  the  internal  wetted  surfaces  of  the  ramjet  become  exposed. 
Subsequently,  when  the  solid  rocket  propellant  is  exhausted,  the  engine  is  operated  in 
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ramjet  mode.  Additional  approaches  are  to  use  other  air  breathing  propulsive  devices 
such  as  turbojets  or  turbofans  in  order  to  achieve  flight  velocities  where  ramjets  become 
effective;  with  ensuing  transition  to  ramjet  propulsion  mode  [9]. 

Ramjets  are  also  susceptible  to  a  phenomenon  known  as  unstart [  .  Unstart  (from 
a  started  condition)  occurs  when  a  normal  shock  wave  or  shock  structure  is  developed  in 
front  of  the  ramjet,  thus  establishing  subsonic  flow  upstream  of  the  ramjet  inlet  face  [I2]. 
Generally,  at  a  given  flight  condition,  a  condition  of  unstart  develops  due  to  the  upstream 
propagation  of  a  normal  shock  or  shock  structures  within  the  ramjet  inlet  until  the  shock 
(or  shocks)  reach  the  converging  portion  of  the  inlet  and  subsequently  entirely  disgorge 
from  the  inlet.  The  resulting  subsonic  flow  entering  the  engine  allows  spillage  (reduced 
mass  capture)  at  the  given  engine  operating  conditions  such  that  a  portion  of  the  air 
associated  with  full  mass  capture  now  flows  around  the  inlet.  This  loss  of  mass  flow  rate 
is  directly  detrimental  to  ramjet  thrust  performance.  In  addition,  there  is  a  large  total 
pressure  loss  through  the  strong  shock  system  upstream  of  the  engine. 

Several  design/operability  strategies  are  used  to  prevent  or  control  unstart.  Two  of 
the  most  common  methods  for  preventing  or  controlling  unstart  include:  1)  the 
incorporation  of  a  dump  combustor  and  2)  bleeding  mass  flow  rate  from  the  ramjet  flow 
path  [13].  Both  of  these  strategies  are  addressed  in  the  present  investigation.  A  dump 
combustor  consists  of  a  sudden  area  increase,  or  step,  after  the  exit  of  the  inlet  and  before 
the  heat  addition.  This  abrupt  increase  in  cross-sectional  area  before  heat  addition  creates 
a  highly  energetic  recirculating  region  that  in  effect  contains  the  normal  shock  and 
prevents  it  from  propagating  upstream  and  unstarting  the  ramjet.  The  dump  combustor 
has  the  added  benefit  of  mixing  the  incoming  fuel  with  the  flow  for  more  complete  and 
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efficient  combustion  in  existing  ramjets  [14].  Bleed  in  a  ramjet  refers  to  a  region,  usually 
within  the  inlet  where  sufficient  pressurization  exists,  where  air  mass  flow  is  extracted 
from  the  main  flow.  Bleed  is  commonly  done  through  a  series  of  traps  or  slots 
incorporated  into  the  design  of  the  wetted  surfaces  of  the  ramjet.  The  normal  shock  may 
propagate  upstream  until  it  encounters  the  region  of  bleed  where,  ideally,  it  is  arrested 
and  kept  from  further  propagating  upstream,  hence  preventing  or  delaying  engine  unstart. 

The  layout  of  this  thesis  consists  of  seven  main  parts.  Section  2,  Ramjet  CFD 
Tools  and  Methodology,  describes  the  CFD  code  utilized  and  the  methods  used  to  derive 
the  ramjet  geometries  used  in  this  investigation.  Section  3,  Ramjet  Throttling  Analytical 
Modeling,  describes  a  method  used  to  analytically  determine  the  maximum  and  minimum 
throttle  heat  addition  rates  for  the  ramjet  configuration.  This  method  involves  the  use  of 
isentropic  flow  relations,  normal  shock  equations,  and  Rayleigh  heat  addition.  Section  4, 
Grid  Convergence  Study,  summarizes  a  trade  study  for  three  different  grids  used  for  a 
baseline  ramjet  configuration  and  determines  whether  or  not  the  grid  used  for  this 
analysis  is  adequately  refined.  Section  5,  Turbulence  Model  Case  Study,  compares  the 
baseline  turbulence  model  used  for  this  analysis  against  another  available  turbulence 
model  in  an  effort  to  characterize  turbulence  model  impact  on  the  unstart  phenomenon. 
Section  6,  Axisymmetric  Inviscid  Wall  Results,  presents  the  results  of  the  inviscid  wall 
analysis  utilizing  the  conventional  (no-dump  combustor)  design  as  well  as  the  dump 
combustor  design.  Section  7,  Axisymmetric  Viscous  Wall  Results,  presents  the  results  of 
the  viscous  walls  analysis  utilizing  only  the  conventional  design.  Section  8,  Summary 
and  Conclusions,  is  a  synopsis  of  the  major  results. 
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2.  RAMJET  CFD  TOOLS  AND  METHODOLOGY 

2.1.  VULCAN  CFD  CODE 

The  VULCAN  (Viscous  Upwind  algorithm  for  Complex  flow  ANalysis)  CFD 
code  version  6.0.2  1151  was  utilized  for  all  the  CFD  results  presented  in  this  work.  In  this 
investigation,  the  viscous  axisymmetric  Navier-Stokes  equations  are  solved  for  both 
inviscid  wall  and  viscous  wall  treatments  of  the  wall  boundaries  in  the  ramjets.  Here 
inviscid  wall  analysis  refers  simply  to  the  use  of  an  inviscid  (slip)  wall  boundary 
condition  for  all  solid  surface  boundaries  of  the  ramjet.  The  viscous  (no-slip)  wall 
analysis  mandates  a  no-slip  boundary  condition  where  the  velocity  at  the  walls  is  defined 
to  be  zero,  with  a  resulting  large  gradient  in  velocity  near  the  walls.  All  solutions  utilized 
adiabatic  walls. 

A  thermally  perfect  gas  mixture  with  76.86%  (by  mass)  nitrogen  and  23.14%  (by 
mass)  oxygen  was  used  to  simulate  the  inflow  air  properties.  Atmospheric  (ambient) 
temperature  and  pressure  corresponded  to  a  standard  altitude  [16]  of  10,000  meters  (223.15 
degrees  Kelvin  and  26,436.3  Pascals,  respectively).  A  free-stream  (flight)  Mach  number 
of  3.0  was  used  for  all  calculations. 

VULCAN  was  run  using  a  time  accurate  distributed  approximating  function 
(DAF)  with  a  dual-time  stepping  scheme.  For  this  scheme,  a  time  step  of  10  nanoseconds 

o 

(1.0x10'  s)  was  used  with  five  sub -iterations  between  the  time  steps.  For  each  sub¬ 
iteration,  a  Courant-Friedrichs-Lewy  number  (CFL)  of  1.0  was  used.  Wall  matching 
functions  were  used  when  the  no-slip  boundary  condition  was  imposed;  grid  refinement 
for  the  configuration  with  viscous  walls  was  not  of  predominant  importance  in  this  study 
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which  largely  focused  on  results  for  inviscid  walls.  The  Menter  Shear  Stress  Transport 
two  equation  (Menter  SST)  model  1171  was  selected  as  the  turbulence  model  for  all 
simulations,  based  on  previous  simulations  and  studies.  Although  the  Menter  SST  model 
is  considered  a  k-co  model,  it  does  blend  the  k-co  (near  the  wall)  with  the  k-s  formulation 
(in  the  free-stream).  A  free-stream  turbulence  intensity  and  turbulent  viscosity  ratio  of 
5%  and  10%  respectively  were  used.  Appendix  A  contains  an  example  VULCAN  input 
deck  used  for  the  viscous  wall  geometry. 

2.2.  GEOMETRY  AND  GRID  DEFINITION 

The  boundary  points  that  defined  the  border  of  the  blocks  that  compose  the 
geometry  of  the  ramjet  were  specified  in  MATLAB  (see  Figure  2.3).  The  boundary 
points  of  the  contoured  wetted  surface  of  the  ramjet  were  calculated  using  cubic  splines 
fitted  between  three  points.  This  method  was  compared  to  results  obtained  with  the 
method  of  characteristics  and  no  significant  difference  was  found  for  this  ramjet 
configuration.  Therefore  the  cubic  spline  method  was  utilized  here  due  to  its  relative 
simplicity.  Figure  2.1  shows  a  comparison  of  the  wall  contours  for  the  diverging  exit 
nozzle  when  using  the  method  of  characteristics  and  cubic  spline  method.  These 
boundary  points  were  then  imported  into  Gridgen  which  allows  the  internal  grid  to  be 
created.  Different  grids  were  constructed  for  both  the  inviscid  and  viscous  wall 
configurations.  For  the  inviscid  wall  internal  grid,  all  of  the  grid  points  were  relatively 
evenly  spaced  throughout  the  domain.  For  the  viscous  wall  grid,  the  internal  grid  points 
were  clustered  at  the  wall  to  capture  the  boundary  layer  effects.  For  this  viscous  wall 
grid,  a  hyperbolic-tangent  distribution  was  used  to  distribute  the  grids  in  the  radial 
direction  starting  at  the  wall  with  a  grid  spacing  of  0.0000025  meters  (2.5pm)  and  then 
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increasing  toward  the  center-line  of  the  ramjet.  The  axial  locations  of  the  grid  points  were 
relatively  evenly  distributed  for  both  the  inviscid  and  viscous  wall  grids. 


Method  of  Characteristics  vs.  Cubic  Spline  Contour 
of  Minimum  Length  Diverging  Nozzle 


Figure  2.1:  Comparison  Between  Boundary  Contours  of  Diverging  Exit  Nozzle  Using 
Method  of  Characteristics  and  Cubic  Spline  Method 


The  area  ratios  used  for  the  inlet  and  exit  nozzles  were  determined  using 
isentropic  quasi-one-dimensional  flow  relations.  The  area  ratio  (local  area  to  throat  area) 
at  a  given  location  is  described  through  the  isentropic  flow  relations,  shown  in  Equation 
(1). 

A  =  M-SKjg 

A*  V  2  / 


M 


(1) 
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Equation  1  can  be  used  to  solve  for  the  Mach  number  at  a  specific  location  given  the 
local  area  ratio.  For  an  inlet  Mach  number  of  3.0,  Equation  (1)  yields  an  area  ratio  of  4.23 
(i.e.  the  inlet  area  is  4.23  times  larger  than  the  throat  area).  However,  per  conventional 
ramjet  design,  it  is  undesirable  to  choke  the  flow  at  the  throat  of  the  inlet  nozzle. 
Therefore  the  area  ratio  for  all  the  geometries  analyzed  has  a  smaller  area  ratio 
(corresponding  to  a  larger  throat  area)  of  1.61.  This  was  determined,  after  CFD 
experimentation,  to  provide  an  adequate  contraction  without  choking  the  flow  in  the  inlet. 
The  inlet  and  exit  radii  for  the  conventional  design  were  set  at  15.625  mm,  making  the 
inlet  area  equal  to  767  mm"  and  the  inlet  throat  area  equal  to  476  mm".  The  dump 
combustor  design  had  the  same  inlet  radius  (15.625  mm),  but  a  larger  exit  radius  of 
23.4375  mm,  due  to  the  area  increase  from  the  dump  combustor,  making  the  exit  area 
equal  to  1,726  mm".  The  exit  nozzle  area  ratio  (Ae/A  )  for  both  ramjet  designs  was  set  at 
1.34,  making  the  exit  throat  areas  equal  to  572  mm"  and  1,288  mm"  for  the  conventional 
and  dump  combustor  designs  respectively.  This  exit  nozzle  area  ratio  yields  two  solutions 
for  Equation  (1);  a  subsonic  solution  with  Mach  number  of  0.50  going  into  the  exit  nozzle 
just  after  combustion  and  a  supersonic  solution  with  Mach  number  of  1.70  at  the  exit  of 
the  ramjet  engine.  This  exit  nozzle  area  ratio  was  used  for  all  of  the  geometries  within 
this  analysis  and  proved  adequate  for  the  purposes  of  the  present  investigation.  With  the 
given  inlet  radius  set  for  both  geometries  (conventional  and  no-dump)  and  the  given 
ambient  conditions  at  the  flight  altitude  (10  km)  the  mass  flow  rate  entering  both  the 
ramjet  designs  was  fixed  at  0.284  kg/s. 

Figure  2.2  shows  the  outline  of  the  final  geometries  analyzed  for  both  the 
conventional  and  dump  combustor  designs. 
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Conventional  Design 


Dunp  Combustor  Design 


Axial  Location  [m] 


Figure  2.2:  Outline  of  Conventional  and  Dump  Combustor  Ramjet  Geometries  Analyzed 


It  can  be  seen  that  the  dump  combustor  geometry  is  0.05  meters  longer  than  the 
conventional  design  (total  length  of  0.25  m).  This  extra  length,  associated  with  the 
constant  area  portion,  allows  the  flow  more  time  to  settle  after  the  dump.  The  exit  lip  of 
both  geometries  has  an  angle  of  approximately  20°  from  the  center-line  direction.  This 
allows  the  pressure  gradient  to  continually  decrease  all  the  way  to  the  exit  plane  of  the 
ramjet. 

2.3.  PARALLELIZATION  STRATEGY 

VULCAN  is  able  to  use  multiple  processors  while  running  a  single  geometry, 
hence  increasing  solution  speed.  The  grid  is  split  into  blocks  and  the  compilation  of  those 
blocks  makes  up  the  overall  geometry.  VULCAN  allocates  one  processor  per  block. 
Figure  2.3  shows  both  the  block  distributions  for  both  conventional  and  dump  combustor 
geometries.  Each  block  has  the  same  number  of  grid  points  for  load  balancing 
(efficiency).  More  blocks  are  used  in  the  diverging  section  of  the  inlet  nozzle  where  the 
grid  is  clustered  (in  order  to  capture  the  large  gradients  associated  with  the  normal  shock 
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that  is  positioned  in  this  component  for  started  flow).  Figure  2.3  also  shows  the  location 
of  the  blocks  where  uniform  heat  addition  (per  volume)  was  specified. 
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Figure  2.3:  Outline  Including  Blocks  and  Heat  Addition  Locations  of  Conventional  and 
Dump  Combustor  Ramjet  Geometries  Analyzed 


2.4.  POST  PROCESSING  RESULTS 

MATLAB  was  used  for  post  processing  all  data  obtained  from  the  VULCAN 
CFD  code.  The  data  collected  from  VULCAN  includes  static  density,  u-velocity  (axial 
component),  v-velocity  (radial  component),  static  pressure,  static  temperature,  Mach 
number,  laminar  viscosity,  eddy  viscosity  ratio,  total  enthalpy  per  mass,  and  the  ratio  of 
specific  heats.  MATLAB  was  then  used  to  calculate  total  entropy  and  entropy  generation 
due  to  friction,  convective  heat  transfer,  ideal  heat  transfer,  non-ideal  heat  transfer,  and 
shocks,  based  on  published  methodologies.  The  MATLAB  code  used  for  post  processing 
can  be  found  in  Appendix  B  along  with  the  details  of  entropy  generation  calculations  in 


Appendix  C. 
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3.  RAMJET  THROTTLING  ANALYTICAL  MODELING 

3.1.  THEORY 

A  quasi-one-dimensional  solver  was  developed  to  predict  the  maximum  and 
minimum  heat  addition  rates  for  the  generic  defined  ramjet.  Minimum  throttle  operation 
was  defined  by  placing  a  normal  shock  within  the  ramjet  at  the  maximum  area 
downstream  of  the  first  throat  (but  upstream  of  heat  addition).  Similarly,  maximum 
throttle  operation  was  defined  by  locating  the  normal  shock  at  the  inlet  throat  (minimum 
area).  In  the  analytical  modeling,  the  flow  is  considered  to  be  isentropic  everywhere 
except  through  shocks  and  within  the  heat  addition  zone. 

Given  the  inlet  conditions  and  the  areas  throughout  the  ramjet,  isentropic  flow 
relations,  normal  shock  equations,  and  the  Rayleigh  flow  1 18  equations  were  then  used  to 
determine  the  amount  of  heat  required  for  a  specified  normal  shock  location  within  the 
diverging  part  of  the  inlet  of  the  ramjet  (i.e.  at  maximum  and  minimum  throttle  operation 
points).  Figure  3.1  shows  the  conventional  geometry  used  in  order  to  calculate  the 
corresponding  heat  addition  rates  associated  with  minimum  and  maximum  throttle 
position,  along  with  the  stations  used  in  the  analysis  involving  the  quasi-one-dimensional 


solver. 


14 


CM 

A 

£ 

p-, _  X. 

6 

r  1 

1''' _ l_l _ l__L. 

1  1  1  '  '  '  ' 

1  1 

1  1  1 

1  1 

'  11  11  . . .  i  i  i  i  . . . i . * 

0  0.05  0.1  0.15  0.2  0.25 

Axial  Location  [m] 


Figure  3.1:  Conventional  Geometry  Used  to  Determine  Maximum  and  Minimum  Throttle 
Heat  Additions  as  Well  as  the  Station  Locations  Used  in  Quasi-One-Dimensional  Solver 


By  starting  calculations  at  the  inlet,  i.e.  station  “0”,  and  progressing  through  the 
ramjet  to  the  exit  plan,  at  station  “6”,  the  heat  addition  rate  corresponding  to  maximum 
throttle  was  calculated  by  the  following  procedure:  First  the  area  ratio  was  calculated 
given  the  Mach  number  at  station  “0”  using  the  isentropic  flow  relation  shown  in 
Equation  (1). 

a,  m 

A*  l  2  /  M„  ^  ' 

A  (the  reference  throat  area)  is  then  calculated,  for  the  given  Ao  (capture  and  inlet  face 
area).  A  new  area  ratio  is  then  calculated  at  station  “1”  and  is  represented  by  Equation 
(2). 

Ai  =  (VHCiFi)  (l+JT^M12)2(Y~i:>  m 

A*  \  2  /  Mi  ^  ' 

Equation  2  is  then  solved  for  Mi  which  represents  the  supersonic  Mach  number  at  station 
“1”,  which  is  the  physical  throat  of  the  inlet.  Using  the  normal  shock  relations  and 
assuming  that  a  normal  shock  stands  at  station  “1”  (as  defined  by  maximum  throttle 
condition);  Equation  (3)  gives  the  Mach  number  just  after  the  shock. 
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Mlsub  represents  the  subsonic  Mach  number  at  station  “1”.  The  total  pressure  is 
conserved  from  the  inlet  face  to  the  normal  shock,  i.e.  p0o  =  p0  .  The  total  pressure  after 
the  normal  shock,  represented  by  p0  is  given  by  Equation  (4). 


Y+l 


Mi 
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ty-l 


Po. 
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=  Po, 
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(4) 


With  Mlsubknown,  Equation  (5)  can  be  solved  to  give  a  new  value  for  the  reference 
throat  area  A  (downstream  of  the  shock)  that  can  be  used  along  with  the  given  area  at 
station  “2”  to  calculate  a  new  area  ratio  at  station  “2”. 


Ai 

A* 


y+l 

2(y— l) 


Y+l 


(v5Aiuk!)2<v_,> 


xsub 


(5) 


This  area  ratio  for  station  “2”  can  in  turn  be  used  to  solve  for  the  Mach  number  at  station 
“2”  (entrance  to  the  combustor)  by  using  Equation  (6). 

Y+l 

— Y+1  ( 2)2(y-1) 

\  2(y-1)  (1+~M2  ) 


M? 


(6) 


The  flow  properties  at  station  “3”  are  the  same  as  station  “2”  given  the  geometry  and 


assumptions  of  isentropic  flow  in  that  region,  i.e.  p0i  =  p0„  =  p0„.  Rayleigh  flow 

isub  z  J 

analysis  [04],  which  is  inviscid  flow  through  a  constant  area  duct,  can  then  be  applied  to 
the  flow  using  Equation  (7)  for  the  static  pressure  change.  Equation  (8)  can  then  be 
solved  for  intermediate  function  f  in  terms  of  the  demanded  total  temperature  change 
through  the  combustor  (determined  by  the  heating  rate). 


P4  =  Ps 


i+ym| 

i+ym| 


(7) 


To. 


(8) 
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where, 


f4  =  (T^j)2(1+1rMJ) 


and 


f=  =  (4k)2(1  +  IrMi> 


u+ym§v 

The  Mach  number  at  station  “4”  can  be  found  by  solving  Equation  (9)  for  the  subsonic 
value  and  assuming  choked  flow  such  that  %  =  — . 

A*  Ag 


y+i 


A4  _  /y+n 

A*  V  2  / 


m4 


(9) 


The  total  pressure  at  station  “4”  can  then  be  calculated  using  Equation  (10). 


Y 

\r- 1 


(10) 


P04  =P4(l  +  ^M|)i 
This  total  pressure  value  is  conserved  through  the  nozzle  and  hence  P04  =  P0g  =  P0g. 
Similarly,  total  temperature  is  conserved  such  that  Tq4  =  T0s  =  T0g.  Equation  (11)  can 
then  be  solved  for  the  supersonic  value  of  Me  again  assuming  choked  flow  at  nozzle 
throat  and  isentropic  flow. 


y+i 


Ag  _  fY+l\ 
A*  V  2  / 


m6 


(11) 


Given  Me,  T0&,  and  P0g,  all  the  exit  plane  quantities  can  be  solved.  Using  the  mass  flow 
rate  entering  the  ramjet,  the  heat  addition  rate  is  related  to  the  change  in  total  temperature 
through  the  combustor  using  Equation  (12). 

Q  =rh(|-l)cPT0o  (12) 

Minimum  throttle  heat  addition  can  be  found  in  a  similar  manner  by  stepping 
through  the  ramjet,  with  the  exception  that  the  normal  shock  is  assumed  to  be  within  the 
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constant  area  portion  of  the  ramjet  between  stations  “2”  and  “3”.  This  simplifies  the 
calculations  because  the  flow  properties  at  the  shock  location  are  the  same  as  at  the  inlet 
face.  This  eliminates  the  need  to  perform  any  calculations  involving  the  inlet.  The  quasi- 
one-dimensional  solver  MATLAB  code  used  for  this  part  of  the  work  is  included  in 
Appendix  D. 

3.2.  ISSUES  INVOLVING  HEAT  ADDITION  MODELING  FOR  ENGINE 
THROTTLING  MANUEVERS 

The  throttling  changes  (maneuvers)  performed  are  referred  to  in  this  work  as  step 
maneuvers  because  of  the  manner  in  which  heat  rate  is  changed  in  order  to  simulate 
throttle  (change)  the  combustion-generated  heat  release.  VULCAN  allows  a  fixed  amount 
of  heat  rate  to  be  specified  within  a  defined  zone  but  does  not  allow  variations  of  that 
heat  addition  rate  within  a  run.  Therefore  the  analysis  modeled  steady-state  at  one  heat 
addition  rate  and  then  changed  this  rate  instantaneously  to  a  new  heat  addition  rate, 
thereby  creating  an  effective  “step”  in  the  amount  of  heat  rate  specified  within  VULCAN. 
While  this  method  of  heat  addition  allows  for  the  most  extreme  transients  to  be  observed 
in  a  conceptual  ‘slam’  throttling  event,  one  must  keep  in  mind  that  the  transients 
observed  in  this  study  are  much  faster  than  practically  observed  in  a  realistic  throttling 
event.  Real  time  lags  in  throttling  (even  ‘slam’  throttling)  are  due  to  a  number  of  factors 
including  the  method  of  fuel  delivery  and  chemical  kinetics  for  the  burn  rates  of  the  fuel. 
However,  it  is  noted  that  developing  technology  may  allow  for  significantly  more  rapid 
throttling  maneuvers  to  take  place  in  future  aircraft  than  is  currently  observed. 
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3.3.  STEADY-STATE  THROTTLING  STUDY 

Steady-state  for  the  entirety  of  this  thesis  is  defined  such  that  the  transients  within 
the  mass  flow  rate  traced  throughout  the  ramjet  are  not  fluctuating  beyond  reasonable 
limits,  the  Mach  number  contours  must  not  be  moving  or  else  fluctuating  with  an 
oscillatory  behavior  within  reasonable  limits,  and  also  by  observing  the  residual  behavior 
across  a  run.  These  conditions  for  steady-state  were  verified  by  observing  the  flow-field 
through  transients,  traced  mass  flow  rate  through  the  ramjet,  and  Mach  number  contours 
every  0.0001  seconds.  These  results  have  been  compiled  into  videos  for  better 
visualization  of  the  transient  behavior  during  throttling  events.  It  should  be  noted  that 
large  recirculation  zones  in  these  types  of  flows  (i.e.  in  the  dump  combustor)  are 
inherently  quasi-steady,  i.e.  generally  exhibit  some  inherent  minor  unsteadiness. 

The  throttling  maneuvers  are  the  same  for  both  ramjet  configurations  and  consist 
of  establishing  a  steady-state  condition  where  no  heat  is  added  (referred  to  in  this  work  as 
“tare”),  then  adding  heat  such  that  a  steady-state  maximum  throttle  condition  is  met 
(referred  to  as  “maximumi”;  the  subscript  one  donating  the  first  maximum  throttle 
setting),  then  subsequently  reducing  the  heat  to  a  corresponding  steady-state  minimum 
throttle  condition  (simply  referred  to  as  “minimum”)  ,  and  then  finally  increasing  the  heat 
addition  such  that  again  a  maximum  throttle  condition  is  met  and  allowed  to  reach 
steady-state  (similarly  referred  to  as  “maximurm”).  Both  the  maximumi  and  maximum2 
throttle  settings  correspond  to  the  same  amount  of  heat  addition  but  are  distinguished  by 
the  subscript  for  ease  of  reference. 

The  reduction  in  heat  rate  from  maximumi  throttle  setting  to  the  minimum  throttle 
setting  will  be  referred  to  as  the  step  down  throttling  maneuver.  Similarly,  changing  the 
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heat  addition  rate  from  the  current  minimum  throttle  setting  to  the  maximum^  throttle 
setting  will  be  referred  to  as  the  step  up  maneuver. 

Using  the  quasi-one-dimensional  solver,  it  was  determined  that  the  nominal 
analytical  minimum  and  maximum  throttle  settings  corresponded  to  approximately  50.4 
and  144.0  kilowatts,  respectively.  These  quantities  were  then  used  as  an  initial  iteration 
for  finding  the  actual  minimum  and  maximum  throttle  positions  based  upon  observations 
of  steady-state  normal  shock  locations  as  found  from  CFD-generated  data. 
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4.  GRID  CONVERGENCE  STUDY 

A  grid  convergence  study  utilizing  inviscid  wall  treatment  was  conducted  in  order 
to  ensure  that  the  results  were  sufficiently  grid-independent  and  also  to  optimize 
computational  efficiency.  Three  grids,  coarse,  medium,  and  fine  were  developed  and 
compared.  The  number  of  blocks  that  composed  the  geometry  remained  the  same  in  all 
three  grids;  the  only  thing  that  changed  was  the  number  of  nodes  per  block.  Table  4.1 
summarizes  the  number  of  nodes  per  block  for  each  grid. 


Table  4.1:  Summary  of  Grid  Resolution  for  the  Grids  Utilized  in  the  Grid  Convergence 

Study 


Grid 

Coarse 

Medium 

Fine 

Nodes  (x,y) 

101,51 

201,101 

301,151 

Blocks 

12 

12 

12 

Total  Nodes 

61,251 

242,501 

543,751 

VULCAN  was  run  with  a  heat  addition  of  144,000  kilowatts,  the  theoretical 
maximum  throttle  setting  from  the  quasi-one-dimensional  solver,  as  a  test  basis  for  all 
three  grids.  Figure  4.1  shows  the  center-line  Mach  number  throughout  the  ramjet  for  all 
three  grids  after  100,000  iterations  or  0.0010  seconds  after  initializing  the  flow  in 


VULCAN. 
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Center  Line  Mach  Number  vs.  Axial  Location 


Figure  4.1:  Center- Line  Mach  Number  Throughout  Ramjet  Utilizing  all  Grids 


It  can  be  seen  that  the  same  trends  are  observed  for  all  three  grids  including  the 
transient  location  of  the  normal  shock  at  approximately  0.125  meters  from  the  inlet.  The 
sudden  drops  in  Mach  number  in  the  inlet  portion  of  the  ramjet  before  the  normal  shock 
correspond  to  weak  oblique  shocks  converging  at  the  center-line. 

Figures  4.2,  4.3,  and  4.4,  show  the  Mach  number  contours  of  a  portion  of  the  grid 
at  the  normal  shock  location  near  the  center-line  of  the  ramjet  corresponding  to  Figure 


4.1. 
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Figure  4.2:  Coarse  Grid  at  Normal  Shock  Location  Near  Center-Line 
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Figure  4.3:  Medium  Grid  at  Normal  Shock  Location  Near  Center-Line 
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Axial  Location  [m] 

Figure  4.4:  Fine  Grid  at  Normal  Shock  Location  Near  Center-Line 


It  can  be  seen  that  while  the  normal  shock  is  approximately  at  the  same  location 
for  each  grid,  the  shock  width  increases  as  the  grid  spacing  increases.  There  is  also  an 
unphysical  jump  in  the  mass  flow  rate  which  can  be  observed  in  Figure  4.5  which  shows 
the  normalized  mass  flow  rates  throughout  the  ramjet  for  all  three  grids  corresponding  to 
Figures  4.2,  4.3,  and  4.4. 

Noting  that  the  spike  in  mass  flow  rate  corresponds  to  the  normal  shock  location 
and  the  magnitude  of  the  spike  decreases  with  increasing  grid  resolution,  the  cause  of  this 
spike  is  most  likely  due  to  grid  resolution  across  the  normal  shock.  It  is  important  to  note 
that  at  this  point  in  the  solution  which  corresponds  to  a  snapshot  at  0.0010  s,  a  steady- 
state  has  not  been  established  and  that  the  normal  shock  is  propagating  upstream  (for  a 
throttle  up  maneuver).  It  was  observed  that  when  the  normal  shock  is  traveling  upstream 
during  the  transients,  mass  flow  rate  tended  to  drop  downstream  of  the  normal  shock 
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location.  The  opposite  of  this  was  observed  as  well,  i.e.  as  the  normal  shock  was  traveling 
downstream  in  a  throttle  down  maneuver,  the  mass  flow  rate  tends  to  increase 
downstream  of  the  normal  shock. 


Figure  4.5:  Normalized  Mass  Flow  Rate  Showing  Normal  Shock  Location,  Spike  in  Mass 
Flow  Rate  at  the  Normal  Shock  Location,  and  Sudden  Drop  in  Mass  Flow  Rate  After  the 

Normal  Shock  for  all  Three  Grids 
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While  all  three  grids  generated  roughly  the  same  Mach  number  contours  and 
provided  the  same  normal  shock  location  at  100,000  iterations,  the  medium  grid  was  used 
for  the  entirety  of  the  analysis  and  was  chosen  as  an  acceptable  trade-off  between 
computational  efficiency  and  accuracy  of  the  results. 
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5.  TURBULENCE  MODEL  CASE  STUDY 

The  objective  of  this  part  of  the  study  was  to  examine  the  impact  of  utilizing  a 
different  turbulence  model  with  emphasis  on  comparing  their  relative  susceptibilities  to 
unstart.  Additional  grid  blocks  were  therefore  added  to  the  configuration  external  to  and 
upstream  of  the  internal  flow  path  in  order  to  allow  unstarted  flow  modeling  (upstream 
and  external  to  the  engine).  Specifically,  the  k-co  model  with  the  Pope  correction,  referred 
to  as  “k-co  Pope”,  is  compared  to  the  base-line  Menter  SST  model  which  is  used  for  the 
bulk  of  the  work  presented  in  this  thesis.  A  major  dilemma  when  modeling  the  ramjet 
with  viscous  wall  treatment,  i.e.  with  the  no-slip  boundary  condition  implemented,  has 
been  the  formation  of  a  large  separated  boundary  layer  upstream  of  the  heat  addition  area. 
This  separated  recirculating  region  is  seen  by  the  main  flow  as  a  flow  blockage,  hence 
effectively  altering  the  geometry  of  the  ramjet.  Given  time,  this  recirculating  region 
propagates  upstream  and  carries  with  it  a  lambda  shock  near  the  walls,  which  transforms 
into  a  normal  shock  toward  the  center-line  of  the  ramjet.  As  soon  as  this  recirculating 
region  reaches  the  throat  of  the  inlet,  incipient  unstart  has  occurred  due  to  the  instability 
of  a  normal  shock  wave  within  a  converging  duct.  An  unstart  condition  is  quickly 
reached  after  this  point  where  a  normal  shock  is  established  upstream  of  the  inlet  of  the 
ramjet.  At  this  point,  subsonic  flow  is  entering  the  engine  and  thus  full  mass  capture  is  no 
longer  possible  (i.e.  the  engine  will  spill  mass).  An  unstarted  engine  produces 
considerably  less  thrust  due  to  the  spillage  and  the  larger  total  pressure  drops  associated 
with  the  stronger  shock  and  considerable  work  has  gone  into  preventing  and  controlling 
unstart.  This  part  of  the  study  therefore  focuses  on  comparing  the  transient  behaviors  of 
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the  flow  in  the  baseline  configuration  as  it  unstarts,  when  using  two  different  turbulence 
models. 

5.1.  MENTER  SST  MODEL 

Due  to  the  unphysical  nature  of  heating  the  flow  (instantaneous  heat  addition) 
very  large  transients  were  observed  during  large  throttle  (heat  rate)  changes.  Figure  5.1 
shows  the  oscillatory  behavior  of  the  normal  shock  before  unstart  occurs. 
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Figure  5.1:  Transient  Mach  Number  Contours  During  Unstart  Event  for  the  Menter  SST 

Model 


Figure  5.2  shows  the  final  Mach  number  contour  after  unstart  at  a  time  of  0.0050 
seconds.  Notice  the  Mach  number  increase  across  the  burner  as  well  as  the  choke 


condition  for  the  exit  nozzle  throat. 
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Figure  5.2:  Mach  Number  Contours  of  Unstart  Condition  at  a  Time  of  0.0050  Seconds 

Utilizing  the  Menter  SST  Model 


Figure  5.3  shows  the  transient  thrust  history  utilizing  the  Menter  SST  model.  The 
moment  of  unstart  can  be  seen  by  the  sudden  decrease  in  thrust  at  approximately  0.0035 
seconds.  The  negative  value  of  thrust  after  unstart  is  evidence  of  the  detrimental  effects 
of  unstart  on  performance. 


Thrust  vs.  Time 


Figure  5.3:  Transient  Thrust  History  for  the  Menter  SST  Model 
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5.2.  k-w  POPE  MODEL 

The  k-G)  Pope  model  was  then  used  for  comparison  in  order  to  see  if  the  transient 
behavior  of  the  unstart  event  was  significantly  different.  Figure  5.4  shows  the  oscillatory 
behavior  of  the  shock  system  for  the  k-co  Pope  model,  which  is  quite  similar  to  Figure 
5.1.  By  observation  of  relative  lambda  shock  locations  between  the  two  models,  it 
appears  that  the  oscillatory  behavior  is  more  severe  for  the  k-co  Pope  model. 
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Figure  5.4:  Transient  Mach  Number  Contours  Before  Unstart  for  the  k-w  Pope  Model 


Figure  5.5  shows  the  final  Mach  number  contours  utilizing  the  k-co  Pope  model  at  a 


time  of  0.0050  seconds  corresponding  to  Figure  5.2.  Again,  unstart  has  been  established. 
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Figure  5.5:  Mach  Number  Contours  of  Unstart  Condition  at  a  Time  of  0.0050  Seconds 

Utilizing  the  k-co  Pope  Model 


Figure  5.6  shows  the  transient  thrust  history  utilizing  the  k-co  Pope  model.  The 
moment  of  unstart  can  again  be  seen  by  the  sudden  decrease  in  thrust  at  approximately 
0.0035  seconds.  The  portions  of  the  graph  that  are  below  -60  Newtons  on  the  y-axis  were 
cut  short  for  scaling  purposes. 


Thrust  vs.  Time 


Figure  5.6:  Transient  Thrust  History  for  the  k-co  Pope  Model 
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Figure  5.7  shows  a  comparison  of  the  mass  flow  rate  through  the  ramjet  after 
unstart  between  the  two  models  (still  within  the  transient  downstream  of  the  normal 
shock).  For  both  models,  nearly  forty  percent  of  the  incoming  mass  is  lost  to  spillage. 
While  the  k-co  Pope  model  has  a  slightly  less  mass  flow  rate  throughout  the  ramjet,  both 
models  show  the  same  trends. 


Time  =  0.0050  s  m  =  0.284  kg/s 


Figure  5.7:  Normalized  Mass  Flow  Rate  Throughout  Ramjet  After  Unstart  for  Both  the 
Menter  SST  and  k-co  Pope  Turbulence  Models 
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The  k-co  Pope  model  was  found  to  be  negligibly  different  from  the  Menter  SST 
model  in  terms  of  transient  behavior  modeling  of  unstart.  All  subsequent  results  shown 
are  therefore  produced  using  the  Menter  SST  model. 
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6.  AXISYMMETRIC  INVISCID  WALL  RESULTS 

Both  geometries  (conventional  and  dump  combustor)  were  modeled  using  the  slip 
boundary  condition  at  the  wall  (here  denoted  as  inviscid  wall).  The  conventional  design, 
i.e.  no-dump  combustor,  was  first  analyzed  in  order  to  determine  maximum  and 
minimum  throttle  settings.  This  was  done  through  trial  and  error  using  the  analytical 
results  discussed  earlier  as  initial  guesses  (with  convergence  assessed  by  examining  the 
flow  until  it  reached  steady-state).  However,  the  ramjet  was  initially  modeled  with  no 
heat  addition  in  the  combustor  in  order  to  establish  an  initial  tare  condition.  Mach 
contours  for  this  tare  case  can  be  seen  in  Figure  6.1. 
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Figure  6.1:  Inviscid  Wall  Steady-State  Mach  Number  Contours  for  No  Heat  Addition  on 

Conventional  Design 


Utilizing  this  tare  flow-field  as  an  initial  condition,  the  theoretical  maximum  heat 
addition  rate,  as  determined  from  the  quasi-one-dimensional  (analytical)  solver,  was  used 
as  a  first  iteration  towards  finding  the  true  maximum  throttle  condition  for  the  CFD. 
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Figure  6.2  shows  the  steady-state  Mach  number  contours  for  this  nominal  (analytical) 
maximum  throttle  setting. 


0  0.05  0.1  0.15  0.2  0.25 

Axial  Location  [m] 


Figure  6.2:  Inviscid  Wall  Steady-State  Mach  Number  Contours  for  a  Theoretical 
Maximum  Heat  Addition  of  144  kW  on  Conventional  Design 


It  can  be  seen  that  maximum  throttle  has  not  yet  been  completely  reached  in  the 
CFD  using  the  analytical  theoretical  value  predicted  by  the  quasi-one-dimensional  model. 
A  throttle  setting  which  was  25%  higher,  i.e.  corresponding  to  180  kW  of  heat  addition, 
was  then  tested  and  was  taken  to  be  the  maximum  throttle,  determined  by  observation  of 
the  normal  shock  location.  Figure  6.3  shows  the  steady-state  Mach  number  contours  of 
this  maximum  throttle  setting. 
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Figure  6.3:  Inviscid  Wall  Steady-State  Mach  Number  Contours  for  Maximum  Throttle 
Setting  of  180  kW  Based  Upon  Normal  Shock  Location  on  Conventional  Design 
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It  can  be  seen  that  at  this  practical  maximum  throttle  position,  the  normal  shock  is 
just  downstream  of  the  inlet  throat  (i.e.  results  correspond  to  a  shock  position  which  is 
slightly  downstream  of  the  throat).  This  is  deemed  acceptable  here  as  it  allows  for  more 
stability.  When  the  shock  is  located  exactly  at  the  inlet  throat,  it  is  neutrally  stable;  any 
perturbation  could  potentially  unstart  the  ramjet.  This  maximum  practical  throttle  setting, 
as  determined  for  the  conventional  ramjet,  was  also  used  for  the  dump  combustor 
geometry.  This  was  done  in  order  to  quantify  the  effects  that  the  dump  combustor  had  on 
Mach  number  position,  stability,  and  entropy  generation.  Multidimensional  effects  most 
likely  account  for  the  25%  increase  in  heat  addition  rate  in  the  actual  CFD-determined 
maximum  throttle  setting,  as  compared  to  the  quasi-one-dimensional  (theoretical)  value. 

From  this  maximum  throttle  setting,  the  heat  rate  was  reduced  to  the  theoretically 
determined  minimum  throttle  amount  of  50.4  kilowatts  (quasi-one-dimensional  analysis) 
as  a  first  guess  for  finding  the  actual  CFD  minimum  throttle  heating  rate.  Figure  6.4 
shows  the  Mach  number  contours  for  this  nominal  minimum  throttle  setting. 
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Figure  6.4:  Inviscid  Wall  Steady-State  Mach  Number  Contours  for  a  Theoretical 
Minimum  Heat  Addition  of  50.4  kW  on  Conventional  Design 
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It  can  be  seen  that  the  normal  shock  is  still  well  within  the  diverging  section  of  the 
inlet  and  thus  true  (multi-dimensional  CFD)  minimum  throttle  operation  has  not  yet  been 
reached.  The  heat  rate  was  reduced  to  720  watts  before  the  minimum  throttle  point  was 
reached,  as  based  upon  normal  shock  location  (recall  the  normal  shock  stands  at  the 
combustor  entrance  for  minimum  throttle).  Figure  6.5  shows  the  Mach  number  contours 
of  the  actual  (CFD-determined)  minimum  throttle  setting. 
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Figure  6.5:  Inviscid  Wall  Steady-State  Mach  Number  Contours  for  Minimum  Throttle 
Setting  of  720  W  Based  Upon  Normal  Shock  Location  on  Conventional  Design 


This  actual  minimum  throttle  setting  is  1.4%  of  the  theoretical  throttle  setting 
calculated  from  the  quasi-one-dimensional  solver.  This  is  again  most  likely  due  to 
multidimensional  affects. 

Figure  6.6  summarizes  the  theoretical  and  CFD  based  maximum/minimum 
throttle  settings  for  the  inviscid  wall  analysis. 
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Theoretical  and  CFD  Based  Maximum/Minimum 
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Figure  6.6:  Summary  of  Theoretical  and  CFD  Based  Maximum/Minimum  Throttle  Heat 

Additions  for  Inviscid  Wall  Analysis 


6.1.  STEADY-STATE  THROTTLING 

Steady-state  throttling  was  determined  as  described  in  section  1;  specifically  by 
observation  of  the  traced  mass  flow  rate  through  the  ramjet  and  Mach  number  contours, 
as  well  as  looking  at  the  L2  residual  per  iteration. 

6.1.1.  Conventional  Design.  Figure  6.6  shows  the  L2  of  the  residual  for  the 
inviscid  wall  conventional  design  analysis.  The  leveling  of  the  residual  corresponds  to  a 
nominal  steady-state  convergence  while  the  sudden  spikes  in  residual  after  steady-state 
has  been  reached  correspond  to  throttle  changes.  This  information  is  also  summarized  in 


Table  6.1. 
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Figure  6.7:  L2  of  the  Residual  vs.  Iteration  for  the  Inviscid  Wall  Conventional  Design 


Table  6.1:  Summary  of  Heat  Addition  and  Iteration  Count  for  Inviscid  Wall  Conventional 

Design  Analysis 


Tare 

Maximum  1 
Throttle 

Minimum 

Throttle 

Maximuni2 

Throttle 

Heat  Addition  [kW] 

0.00 

180.00 

0.72 

180.00 

Current  Iterations 

100,000 

500,000 

700,000 

500,000 

Total  Iteration  Count 

100,000 

600,000 

1,300,000 

1,800,000 
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6.1. 1.1  Fluid  dynamics.  Figure  6.8  summarizes  the  steady-state  Mach 
number  contours  achieved  for  the  inviscid  wall  analysis  for  the  conventional  design. 
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Figure  6.8:  Inviscid  Wall  Summary  of  Mach  Number  Contours  at  Major  Throttle  Settings 

for  Conventional  Design 


6.1. 1.2  Performance  and  entropy  results  and  analysis.  The  thrust  data 
achieved  for  the  steady-state  analysis  of  the  inviscid  wall  conventional  design  is 
presented  in  Table  6.2. 


Table  6.2:  Summary  of  Thrust  Values  for  Steady-State  Inviscid  Wall  Conventional 


Design 


Tare 

Maximumi 

Throttle 

Minimum 

Throttle 

Maximurm 

Throttle 

Thrust  [N] 

-6.9 

54.7 

-34.7 

55.0 
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The  thrust  for  the  minimum  throttle  is  less  than  the  tare  case.  This  is  most  likely 
due  to  the  losses  from  the  strong  normal  shock  within  the  flow  path  for  the  minimum 
throttle  setting. 

Figures  6.9,  6.10,  6.11  and  6.12  are  a  breakdown  of  the  total  entropy  and  the 
entropy  generation  due  to  different  loss  mechanisms  (friction,  convective  heat  transfer, 
non-ideal  heat  transfer,  ideal  heat  transfer,  and  shocks)  through  the  ramjet  for  the  tare, 
maximumi,  minimum,  and  maximum2  throttle  settings  respectively.  The  entropy  values 
have  been  scaled  from  the  inlet  entropy  such  that  the  inlet  entropy  value  is  zero.  For  all 
throttle  settings,  entropy  generated  due  to  friction  and  convective  heat  transfer  is 
negligible,  a  direct  consequence  of  the  slip  wall  boundary  condition  and  the  largely  one¬ 
dimensional  nature  of  the  flow.  It  is  interesting  to  note  that  entropy  generated  due  to  non¬ 
ideal  heat  transfer  is  also  very  small  relative  to  ideal  heat  transfer,  implying  that  the 
majority  of  the  heat  transfer  takes  place  ideally.  This  is  not  surprising,  given  the  low 
Mach  number  in  the  heat  release  zone  (where  the  static  temperature  is  not  significantly 
lower  than  the  total  temperature). 

It  is  important  to  note  that  the  overall  entropy  should  never  decrease  throughout 
the  ramjet  for  steady  flow.  Any  such  decrease  in  steady  flow  would  be  a  violation  of  the 
second  law  of  thermodynamics.  Note,  however,  that  entropy  generation  rate  and  hence 
entropy  generation  per  mass  can  and  does  change  within  the  transient  process  in  a 
throttling  event  (i.e.  not  in  a  steady-state  situation).  All  flow  rates,  including  the  mass 
flow  rate  change  markedly,  displaying  oscillatory  behavior,  during  the  transient  event 
itself.  This  sort  of  behavior  will  be  evident  in  some  subsequent  plots  and  results  which 
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are  representative  of  transient  flow-fields.  For  some  cases,  inherent  unsteadiness  in  the 
flow-field  made  an  absolute  steady-state  solution  impractical,  if  not  impossible. 


Entropy  vs.  Axial  Location 


Figure  6.9:  Entropy  Breakdown  for  Inviscid  Wall  Conventional  Design  Tare  Case 


An  example  of  this  (very  slight)  decrease  in  entropy  can  be  seen  in  Figure  6.9 
above,  which  illustrates  the  entropy  breakdown  for  the  tare  case.  However,  note  that  the 
overall  magnitude  of  the  entropy  throughout  the  ramjet  for  the  tare  case  is  extremely 
small  and  coupled  with  the  transient  effects  is  considered  negligible,  especially  when 
compared  to  the  (subsequent)  entropy  distributions  associated  with  the  throttled  cases. 
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Entropy  vs.  Axial  Location 


Axial  Location  [m] 

Figure  6.10:  Entropy  Breakdown  for  Inviscid  Wall  Conventional  Design  Maximum] 

Throttle  Case 


Figure  6.10  provides  the  entropy  distribution  in  terms  of  loss  mechanism  for  the 
inviscid  wall  conventional  design  (maximum]  throttle)  at  a  nominal  steady-state.  This 
figure  clearly  shows  the  unphysical  spike  in  entropy  corresponding  to  the  normal  shock 
location  within  the  diverging  portion  of  the  inlet.  This  spike  in  entropy  is  due  to  the  grid 
resolution  across  the  normal  shock  and  is  associated  with  a  numerical/grid  induced  spike 
in  mass  flow  rate  within  the  CFD  solutions  as  discussed  earlier.  This  spike  in  entropy 
(and  mass)  is  observed  in  all  of  the  results  where  a  normal  shock  is  present  and  will  not 
be  further  discussed  in  subsequent  results  presented  here. 

Figure  6.10  also  shows  the  dominance  of  the  shock-induced  entropy  rise  (first 
rise)  and  the  heat  addition  entropy  rise  (second  rise)  in  this  ramjet  flow.  The  influence  of 
both  friction  and  heat  transfer  are  very  small  across  the  length  of  the  ramjet.  It  can  also  be 
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seen  that  the  total  entropy  and  entropy  from  shocks  is  very  slightly  decreasing  from  the 
normal  shock  location  to  the  exit  of  the  ramjet  (excluding  the  increase  in  total  entropy 
across  the  burner).  Recall  that  while  the  normal  shock  is  traveling  upstream  during  a 
transient  maneuver,  the  mass  flow  (and  hence  other  variables  depending  on  mass  flow) 
decreases  after  the  normal  shock.  This  causes  the  slight  decrease  in  entropy  observed 
after  the  normal  shock,  even  at  this  nominally  steady-state  point.  This  trend  is  observed  in 
almost  all  of  the  following  steady-state  analysis  (the  opposite  effect  will  be  observed 
while  the  normal  shock  is  moving  downstream,  i.e.  an  increase  in  mass  flow  rate  and 
other  variables  is  observed  after  the  normal  shock).  This  is  a  consequence  of  the  limited 
time  and  resources  available  that  would  be  (conceptually)  required  to  reach  a  completely 
steady-state  solution. 

Figure  6.11  shows  the  entropy  generation  history  for  the  minimum  throttle 
position  for  nominal  steady-state.  Here  the  normal  shock  has  moved  downstream  and  is 
contained  within  the  constant  area  portion  of  the  ramjet.  Even  though  the  shock  is 
stronger  and  the  entropy  rise  associated  with  the  shock  is  larger  than  the  maximum 
throttle  cases,  the  overall  entropy  is  significantly  less  when  compared  to  the  maximum 
throttle  cases.  This  is  due  to  the  considerable  decrease  in  heat  addition  rate  and  its 
corresponding  entropy  rise. 

Figure  6.12  is,  as  expected  since  they  correspond  to  the  same  throttle  setting, 
noticeably  similar  to  Figure  6.10  and  is  shown  here  for  completeness. 
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Figure  6.1 1:  Entropy  Breakdown  for  Inviscid  Wall  Conventional  Design  Minimum 

Throttle  Case 
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Figure  6.12:  Entropy  Breakdown  for  Inviscid  Wall  Conventional  Design  Maximuirb 

Throttle  Case 


45 


Figures  6.13,  6.14,  6.15,  and  6.16  show  the  percentage  breakdown  of  entropy  as  a 
percentage  of  the  total  cumulative  entropy  at  the  exit  plane  of  the  ramjet  for  the  tare, 
maximumi,  minimum,  and  maximum2  throttle  settings  respectively  on  the  inviscid  wall 
conventional  design.  In  essence,  these  figures  are  a  reiteration  of  the  exit  plane  values  of 
entropy  corresponding  to  Figures  6.9,  6.10,  6.11,  and  6.12.  Without  heat  addition  (i.e.  in 
the  tare  case),  entropy  generation  due  to  friction,  non-ideal  and  ideal  heat  transfer  is 
virtually  zero. 


Inviscid  Wall  Entropy  Allocation,  Tare 


Figure  6.13:  Percentage  Breakdown  of  Entropy  Due  to  its  Mechanisms  for  Tare  Case  at 
the  Exit  Plane  on  Inviscid  Wall  Conventional  Design 
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For  the  tare  case  (Figure  6.13),  the  only  significant  contributor  to  entropy 
generation  is  due  to  the  relatively  weak  oblique  shocks  within  the  ramjet.  These  shocks 
account  for  97.7%  of  the  total  entropy  generated,  with  entropy  due  to  friction  and 
convective  heat  transfer  accounting  for  the  remaining  2.3%.  However,  it  should  be  noted 
that,  by  comparison,  the  tare  case  only  produces  1.4%  of  the  total  amount  of  entropy 
generated  in  the  maximum2  throttle  case. 


Inviscid  Wall  Entropy  Allocation,  Maximumj 
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Figure  6.14:  Percentage  Breakdown  of  Entropy  Due  to  its  Mechanisms  for  Maximumi 
Throttle  Case  at  the  Exit  Plane  on  Inviscid  Wall  Conventional  Design 
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For  the  maximumi  throttle  case  (Figure  6.14),  it  is  seen  that  80.2%  of  the  entropy 
generation  is  due  to  ideal  heat  transfer  associated  with  the  heat  addition,  while  shocks 
account  for  16.2%,  with  the  other  mechanisms  comprising  the  remaining  3.7%. 


In  viscid  Wall  Entropy  Allocation,  Minimum 


Figure  6.15:  Percentage  Breakdown  of  Entropy  Due  to  its  Mechanisms  for  Minimum 
Throttle  Case  at  the  Exit  Plane  on  Inviscid  Wall  Conventional  Design 


For  the  minimum  throttle  case  (Figure  6.15),  similarly  to  the  tare  case,  shocks  are 
the  largest  contributor  of  the  entropy  generation  (76.8%),  followed  by  convective  heat 
transfer  (21.9%).  By  comparison,  the  minimum  throttle  case  only  generates  41.3%  of  the 
total  entropy  generated  in  the  maximum2  throttle  case. 
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Inviscid  Wall  Entropy  Allocation,  Maximum2 
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Figure  6.16:  Percentage  Breakdown  of  Entropy  Due  to  its  Mechanisms  for  Maximum? 

Throttle  Case  at  the  Exit  Plane  on  Inviscid  Wall  Conventional  Design 

The  maxi  mum  2  throttle  case  (Figure  6.16),  as  expected,  is  very  similar  to  the 
maximumi  throttle  case,  (shown  in  Figure  6.14),  but  is  shown  here  for  completeness.  A 
summary  of  the  entropy  allocation  at  the  exit  plane  of  the  ramjet  for  all  the  throttling 
cases  can  be  found  in  Table  6.3  for  the  inviscid  wall  conventional  design. 
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Table  6.3:  Summary  of  Entropy  Allocation  at  Exit  Plane  of  Ramjet  for  All  Steady-State 


Throttling  Cases : 

for  Inviscid  Wall  Conventional  Design 

Mechanism  [J/kg-K] 

Tare 

Maximum] 

Minimum 

Maximurm 

Friction 

0.01 

0.06 

0.20 

0.08 

Convective  HT 

0.07 

3.83 

23.42 

4.85 

Non-ideal  HT 

0.00 

5.50 

0.05 

5.48 

Ideal  HT 

0.00 

204.00 

1.17 

204.24 

Shocks 

3.48 

41.11 

82.14 

44.22 

Total 

3.56 

254.49 

106.99 

258.87 

6.1.2.  Dump  Combustor.  Figure  6.17  shows  the  L2  of  the  residual  for  the 
inviscid  wall  dump  combustor  design.  This  information  is  also  summarized  in  Table  6.4. 
By  comparison  with  the  inviscid  wall  conventional  design  residual  history  shown  in 
Figure  6.7,  the  dump  combustor  exhibits  much  more  variability  in  the  residual.  This  is 


expected  due  to  the  large  recirculating  region  associated  with  the  dump. 
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Figure  6.17:  L2  of  the  Residual  vs.  Iteration  for  Inviscid  Wall  Dump  Combustor  Design 


Table  6.4:  Summary  of  Fleat  Addition  and  Iteration  Count  for  Inviscid  Wall  Dump 

Combustor  Design 


Tare 

Maximum 

Throttle 

Minimum 

Throttle 

Maximum 

Throttle 

Heat  Addition  [kW] 

0.00 

180.00 

0.72 

180.00 

Current  Iterations 

100,000 

800,000 

500,000 

500,000 

Total  Iteration  Count 

100,000 

900,000 

1,400,000 

1,900,000 

51 


6.1.2.1  Fluid  dynamics.  Figure  6.18  summarizes  the  steady-state  Mach 
number  contours  achieved  for  the  inviscid  wall  analysis  of  the  dump  combustor  design.  It 
can  be  seen  that  at  minimum  throttle,  the  normal  shock  is  completely  swept  out  of  the 
ramjet  exit,  leaving  supersonic  flow  throughout  the  ramjet. 


E 
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Figure  6.18:  Summary  of  Mach  Number  Contours  for  Inviscid  Wall  Dump  Combustor 

Design 


The  sudden  area  change  at  the  dump  location  also  produces  a  noticeable 
recirculation  region  that  is  confined  to  the  region  just  aft  of  the  dump  in  the  tare  and 
minimum  throttle  setting,  but  is  much  more  pronounced  and  extends  into  the  heat 
addition  region  in  the  maximum  throttle  settings  as  shown  in  Figure  6.19. 
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Figure  6.19:  Recirculating  Regions  on  Minimum  and  Maximuni2  Throttle  Settings  for 

Inviscid  Wall  Dump  Combustor  Design 


6.1.2.2  Performance  and  entropy  results  and  analysis.  The  thrust  data 
achieved  for  the  steady-state  analysis  is  contained  in  Table  6.5. 


Table  6.5:  Summary  of  Thrust  Values  for  Steady-State  Inviscid  Wall  Dump  Combustor 


Design 


Tare 

Maximum  i 
Throttle 

Minimum 

Throttle 

Maximurm 

Throttle 

Thrust  [N] 

7.9 

21.4 

7.9 

23.1 

For  the  dump  combustor  case,  the  minimum  throttle  thrust  is  approximately  equal 
to  the  tare  thrust  due  to  the  flow  similarities  between  the  two. 

Figures  6.20,  6.21,  6.22  and  6.23  provide  breakdowns  of  the  entropy  generation 
due  to  different  loss  mechanisms  for  the  tare,  maximum i,  minimum,  and  maximum2 
throttle  settings,  respectively.  While  entropy  generation  due  to  friction  and  convective 
heat  transfer  are  still  small,  the  magnitude  of  the  overall  entropy  generation  is  much 
larger  than  observed  in  the  conventional  inviscid  wall  analysis.  It  is  observed  that  the 
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dump  combustor  has  a  general  trend  to  increase  the  contributions  to  overall  entropy  from 
mechanisms  other  than  shocks  and  heat  transfer  when  compared  to  the  conventional 
inviscid  wall  design. 


Entropy  vs.  Axial  Location 


Axial  Location  [m] 

Figure  6.20:  Entropy  Breakdown  for  Tare  Case  on  Inviscid  Wall  Dump  Combustor 

Design 


Figure  6.20  shows,  for  the  tare  case  (dump  combustor),  a  large  increase  in  entropy 
due  to  shocks  downstream  of  the  dump  and  also  a  slight  increase  in  entropy  due  to 
convective  heat  transfer  (comprising  the  total  heat  transfer  for  the  tare  case).  It  is 
important  to  note  that  for  the  inviscid  wall  dump  combustor  analysis,  steady-state  was 
never  completely  achieved,  again,  due  to  the  limited  time  and  resources  and  the  inherent 
unsteadiness  associated  with  the  large  recirculating  region  in  the  dump  combustor. 
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Figure  6.21:  Entropy  Breakdown  for  Maximumi  Throttle  Case  on  Inviscid  Wall  Dump 

Combustor  Design 


Figure  6.21  shows  that  for  the  maximumi  throttle  case  (dump  combustor)  the 
major  contributor  to  entropy  generation  is  associated  with  the  ideal  heat  transfer  loss 
mechanism.  The  combined  effect  of  friction  and  convective  heat  transfer  for  this  case  is 
over  10  times  greater  than  the  corresponding  effect  in  the  inviscid  wall  conventional 
design  analysis.  Also  note  that  entropy  generated  due  to  non-ideal  heat  transfer  has 
significantly  increased  when  compared  to  the  inviscid  wall  conventional  design  analysis. 
This  can  be  explained  by  the  expansion  of  the  flow  across  the  dump  combustor  and  the 
increase  in  Mach  number  at  which  heat  is  added  (the  higher  Mach  number  causes  higher 


non-ideal  heat  transfer). 
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Entropy  vs.  Axial  Location 
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Figure  6.22:  Entropy  Breakdown  for  Minimum  Throttle  Case  on  Inviscid  Wall  Dump 

Combustor  Design 


Figure  6.22  (minimum  throttle  dump  combustor)  is,  as  expected,  similar  to  Figure 
6.20  for  the  tare  case.  This  is  expected  since  the  flow-fields  are  nearly  identical  in 
features.  The  major  difference  is  the  contribution  to  entropy  generation  from  non-ideal 
and  ideal  heat  transfer  in  the  burner. 

Figure  6.23  (entropy  allocation  for  maximurm  throttle  setting)  is,  again,  similar  to 
Figure  6.21  (entropy  allocation  for  maximum i  throttle  setting). 
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Entropy  vs.  Axial  Location 


Axial  Location  [m] 

Figure  6.23:  Entropy  Breakdown  for  Maximurm  Throttle  Case  on  Inviscid  Wall  Dump 

Combustor  Design 


Figures  6.24,  6.25,  6.26,  and  6.27  show  the  percentage  breakdown  of  entropy  as  a 
percentage  of  the  total  cumulative  entropy  at  the  exit  plane  of  the  ramjet  for  the  tare, 
maximumi,  minimum,  and  maximumi  throttle  settings  respectively  for  the  inviscid  wall 
dump  combustor  design.  In  essence,  these  figures  are  a  reiteration  of  the  exit  plane  values 
of  entropy  corresponding  to  Figures  6.20,  6.21,  6.22,  and  6.23. 
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Figure  6.24:  Percentage  Breakdown  of  Entropy  Generation  Due  to  its  Mechanisms  for 
Tare  Case  on  Inviscid  Wall  Dump  Combustor  Design 


While  the  total  entropy  for  the  inviscid  wall  dump  combustor  tare  case  as  shown 
in  Figure  6.24  is  small  compared  to  the  maximum  throttle  cases,  it  is  approximately  17 
times  larger  than  the  corresponding  tare  case  for  the  inviscid  wall  conventional  design 
analysis.  Again,  in  the  tare  case,  the  only  major  contributor  to  entropy  is  the  relatively 
weak  oblique  shocks  within  the  ramjet.  These  shocks  account  for  84.0%  of  the  total 
entropy  generated,  with  entropy  generated  due  to  friction  and  convective  heat  transfer 
accounting  for  the  remaining  16.0%.  Interestingly,  by  comparison,  the  tare  case  now 
produces  18.9%  of  the  total  amount  of  entropy  generated  in  the  maxi  mum  2  throttle  case, 
compared  to  the  1 .4%  in  the  conventional  inviscid  wall  tare  case. 
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Figure  6.25:  Percentage  Breakdown  of  Entropy  Generation  Due  to  its  Mechanisms  for 
Maximumi  Throttle  Case  on  Inviscid  Wall  Dump  Combustor  Design 


For  the  maximumi  throttle  case  shown  in  Figure  6.25,  it  can  be  seen  that  63.8%  of 
the  entropy  generation  is  due  to  ideal  heat  transfer.  The  entropy  due  to  shocks  is  negative 
due  to  the  nature  of  the  step  up  maneuver  and  also  due  to  unsteadiness  (transient  effects). 
The  combined  effects  of  friction  and  convective  heat  transfer  for  the  inviscid  wall  dump 
combustor  now  accounts  for  14.7%  of  the  total  entropy  generated  within  the  ramjet, 
compared  to  only  1.5%  for  the  corresponding  inviscid  wall  conventional  design  analysis. 
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Figure  6.26:  Percentage  Breakdown  of  Entropy  Generation  Due  to  its  Mechanisms  for 
Minimum  Throttle  Case  on  Inviscid  Wall  Dump  Combustor  Design 

The  inviscid  wall  dump  combustor  minimum  throttle  case  (Figure  6.26)  yielded 
similar  results  to  the  conventional  inviscid  wall  case  with  the  exception  of  non-ideal  heat 
transfer.  Non-ideal  heat  transfer  now  has  a  greater  impact  on  entropy  compared  to  ideal 
heat  transfer.  For  the  minimum  throttle  case,  the  entropy  generated  is  20.3%  of  the  total 
entropy  generated  in  the  maximum2  throttle  case.  This  is  comparable  to  the  tare  case 
(18.9%)  as  expected  due  to  the  absence  of  a  normal  shock  present  in  both  cases. 
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Figure  6.27 :  Percentage  Breakdown  of  Entropy  Generation  Due  to  its  Mechanisms  for 
Maximumo  Throttle  Case  on  Inviscid  Wall  Dump  Combustor  Design 


The  maximumo  throttle  case  (Figure  6.27)  is  very  similar  to  the  maximum i 
throttle  case  (Figure  6.25).  A  summary  of  the  entropy  allocations  at  the  exit  plane  of  the 
ramjet  for  all  the  throttling  cases  can  be  found  in  Table  6.6  for  the  inviscid  wall  dump 


combustor  design. 
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Table  6.6:  Summary  of  Entropy  Allocation  at  Exit  Plane  of  Ramjet  for  All  Steady-State 
Throttling  Cases  for  Inviscid  Wall  Dump  Combustor  Design 


Mechanism  [J/kg-K] 

Tare 

Maximumi 

Minimum 

Maximurm 

Friction 

1.89 

31.74 

2.36 

30.78 

Convective  HT 

7.63 

14.33 

7.41 

19.44 

Nonideal  HT 

0.00 

49.60 

3.48 

47.55 

Ideal  HT 

0.00 

199.85 

1.15 

202.33 

Shocks 

50.12 

-17.51 

49.35 

14.67 

Total 

59.64 

278.00 

63.76 

314.77 

6.2.  TRANSIENT  THROTTLING 

This  section  provides  results  of  transient  throttling  events  for  both  the 
conventional  (no-dump)  and  the  dump  combustor  configurations  with  inviscid  walls. 

6.2.1.  Conventional  Design  This  section  presents  the  fluid  dynamic, 
performance,  and  entropic  analysis  for  the  transient  throttling  maneuvers  performed  on 
the  inviscid  wall  conventional  design. 

6.2.1. 1  Fluid  dynamics.  The  transients  of  the  throttle  step  down  maneuver 

(maximum  to  minimum  throttle)  can  be  seen  in  the  snap  shots  of  the  Mach  number 
contours  captured  in  Figure  6.28.  Each  frame  represents  a  time  step  of  0.0010  seconds  for 
a  total  elapsed  time  of  0.0070  seconds  and  approximately  covers  the  duration  of  the 


transient  maneuver. 
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Figure  6.28:  Mach  Number  Contours  of  Step  Down  Maneuver  for  Inviscid  Wall 

Conventional  Design 


Figure  6.29  shows  the  mass  flow  rate  through  the  ramjet  corresponding  to  Figure 

6.28.  Notice  that  the  normal  shock  locations  in  Figure  6.28  correspond  to  those  in  Figure 

6.29.  Evidence  of  the  increase  in  mass  flow  rate  after  the  normal  shock  during  a  step 
down  maneuver  (normal  shock  propagating  downstream)  is  apparent.  The  magnitude  of 
the  increase  in  mass  flow  diminishes  as  a  steady- state  solution  is  almost  asymptotically 


approached. 
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m  =  0.284  kg/s 


Figure  6.29:  Mass  Flow  Rate  Through  Ramjet  During  Transient  Step  Down  Maneuver 

for  Inviscid  Wall  Conventional  Design 


The  transients  of  the  step  up  maneuver  can  be  seen  in  the  snap  shots  of  the  Mach 
number  contours  captured  in  Figure  6.30.  Each  frame  represents  a  time  step  of  0.0010 
seconds  for  a  total  elapsed  time  of  0.0050  seconds  approximately  covering  the  duration  of 
the  transient  event  similar  to  the  step  down  maneuver.  A  steady-state  solution  seems  to  be 
reached  soon  relative  to  the  step  down  maneuver. 
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Figure  6.30:  Mach  Number  Contours  of  Step  Up  Maneuver  for  Inviscid  Wall 

Conventional  Design 


Figure  6.31  shows  the  mass  flow  rate  through  the  ramjet  corresponding  to  Figure 
6.30.  Notice  again  that  the  normal  shock  locations  in  Figure  6.30  correspond  to  those  in 
Figure  6.31.  Evidence  of  the  decrease  in  mass  flow  rate  after  the  normal  shock  during  a 
step  up  maneuver  (normal  shock  propagating  upstream)  is  apparent.  The  magnitude  of 
the  increase  in  mass  flow  after  the  normal  shock  diminishes  as  a  steady-state  solution  is 
almost  asymptotically  approached  similar  to  the  step  down  maneuver. 
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Figure  6.31:  Mass  Flow  Rate  Through  Ramjet  During  Transient  Step  Up  Maneuver  for 

Inviscid  Wall  Conventional  Design 


6.2.1.2  Performance  and  entropy  results  and  analysis.  Figure  6.32  shows  the 
transient  thrust  behavior  for  the  inviscid  wall  conventional  design  analysis.  Tare  is  first 
specified  and  after  0.001  seconds  has  established  an  approximate  steady-state  flow-field 
(sees  Table  6.2  for  the  exact  steady-state  values  of  thrust  for  all  major  throttle  settings). 
Maximum  heat  addition  is  then  added  and  a  severe  decrease  in  thrust  is  noticed  before  the 


thrust  climbs  and  approaches  a  steady  value  at  0.006  seconds  into  the  analysis.  At  this 
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point  in  time,  the  step  down  maneuver  is  performed  and  a  sudden  slight  increase  in  thrust 
is  experienced  before  the  thrust  quickly  decreases  to  a  steady  negative  value  after  0.013 
seconds.  After  this  steady  minimum  throttle  is  established,  the  step  up  maneuver  is 
performed  at  0.013  seconds  and  approaches  a  steady  value,  similar  to  the  initial 
maximum  throttle,  after  0.018  seconds  and  the  completion  of  the  analysis. 


Thrust  vs.  Time 


Figure  6.32:  Transient  Thrust  Behavior  for  Inviscid  Wall  Conventional  Design 


Figure  6.33  shows  the  transient  behavior  of  the  total  entropy  per  mass  at  the  exit 
plane  of  the  ramjet.  The  trend  for  the  entropy  at  the  onset  of  the  transient  maneuvers  is 
approximately  reversed  from  what  was  observed  in  the  thrust  history  for  the  same 
transient.  Specifically,  when  maximum  throttle  is  initiated  (from  the  tare  flow),  the  thrust 
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value  initially  decreases  and  then  rises.  In  the  case  of  the  entropy,  there  is  a  drastic  rise, 
then  a  rapid  decrease  below  the  initial  value,  followed  by  a  final  rise  to  an  intermediate 
value.  This  opposing  transient  behavior  (between  thrust  and  entropy)  is  also  observed  for 
the  step  down  maneuver;  the  thrust  first  slightly  increased  before  decreasing  while 
entropy  first  decreased,  then  increased.  This  trend  is  expected  in  transient  situations  and 
is  mainly  due  to  mass  flow  rate  issues,  since  thrust  and  entropy  per  mass  respond 
differently  to  mass  flow  rate  transients.  In  addition,  irreversibilities  decrease 
performance. 


Entropy  vs.  Time 


Figure  6.33:  Transient  Total  Entropy  Behavior  for  Inviscid  Wall  Conventional  Design 
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6.2.2.  Dump  Combustor.  This  section  presents  the  fluid  dynamic, 
performance,  and  entropic  analysis  for  the  transient  throttling  maneuvers  performed  on 
the  inviscid  wall  dump  combustor  design. 

6.2.2. 1  Fluid  dynamics.  The  transients  of  the  step  down  maneuver  for  the 
dump  combustor  (inviscid  walls)  can  be  seen  in  the  snap  shots  of  the  Mach  number 
contours  captured  in  Figure  6.34.  Each  frame  represents  a  time  step  of  0.0010  seconds  for 
a  total  elapsed  time  of  0.0050  seconds  approximately  covering  the  duration  of  the 
transient  maneuver.  As  noted  earlier,  the  normal  shock  system  passes  through  the  exit 
nozzle  and  is  no  longer  present  at  the  minimum  throttle. 
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Figure  6.34:  Mach  Number  Contours  of  Step  Down  Maneuver  for  Inviscid  Wall  Dump 

Combustor  Design 
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Figure  6.35  shows  the  mass  flow  rate  through  the  ramjet  corresponding  to  Figure 
6.34.  Evidence  of  the  increase  in  mass  flow  rate  after  the  normal  shock  is  apparent  while 
the  normal  shock  is  still  within  the  ramjet  before  it  passes  through  the  exit  nozzle.  A 
steady-state  flow-field  is  quickly  achieved  in  the  absence  of  the  normal  shock  system. 


m  =  0.284  kg/s 


Figure  6.35:  Mass  Flow  Rate  Through  Ramjet  During  Transient  Step  Down  Maneuver 
for  Inviscid  Wall  Dump  Combustor  Design 
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The  transients  of  the  step  up  maneuver  can  be  seen  in  the  snap  shots  of  the  Mach 
number  contours  captured  in  Figure  6.36.  Each  frame  represents  a  time  step  of  0.0010 
seconds  for  a  total  elapsed  time  of  0.0050  seconds  approximately  covering  the  duration  of 
the  transient  maneuver  similar  to  the  step  down  maneuver.  A  quasi-steady  solution  seems 
to  be  reached  fairly  soon  relative  to  the  step  down  maneuver  which  was  similar  to  the 
results  for  the  inviscid  wall  conventional  design  analysis. 
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Figure  6.36:  Mach  Number  Contours  of  Step  Up  Maneuver  for  Inviscid  Wall  Dump 

Combustor  Design 


Figure  6.37  shows  the  mass  flow  rate  through  the  ramjet  for  the  step  up  maneuver 
corresponding  to  Figure  6.36.  Initially,  drastic  oscillations  are  observed  that  quickly 
diminish  to  reach  a  relatively  steady-state  solution.  Because  there  is  no  distinct 
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propagation  of  a  normal  shock  upstream  of  the  dump  combustor,  the  pronounced 
decrease  in  mass  flow  rate  that  would  normally  be  observed  is  absent. 


m  =  0.284  kg/s 


Figure  6.37:  Mass  Flow  Rate  Through  Ramjet  During  Transient  Step  Up  Maneuver  for 

Inviscid  Wall  Dump  Combustor  Design 
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6.2.2.2  Performance  and  entropy  results  and  analysis.  Figure  6.38  shows  the 
transient  thrust  behavior  for  the  inviscid  wall  dump  combustor  design  analysis.  Tare  is 
first  specified  and  after  0.001  seconds  has  established  a  steady  value  of  thrust  (sees  Table 
6.4  for  the  exact  steady-state  thrust  values  for  all  major  throttle  settings).  Maximum  heat 
addition  rate  is  then  added  and  a  severe  decrease  in  thrust  is  noticed  before  the  thrust 
climbs  and  approaches  a  quasi  steady  value  at  approximately  0.006  seconds  into  the 
analysis.  This  quasi-steady  thrust  value  then  decreases  (with  specified  heat  rate  remaining 
constant)  and  reaches  a  different  lower  steady  value  after  0.009  seconds.  This  behavior  is 
most  likely  due  to  the  inherent  unsteadiness  of  the  dump  combustor  and  the  passing  of  a 
recirculation  region  through  the  ramjet.  At  this  point  in  time,  the  step  down  maneuver  is 
performed  and  a  sudden  slight  increase  in  thrust  is  experienced  before  the  thrust  quickly 
decreases  while  exhibiting  large  oscillatory  behavior  before  reaching  a  steady  value  after 
0.014  seconds  into  the  analysis.  After  this  steady  minimum  throttle  condition  is 
established,  the  step  up  maneuver  is  performed  at  0.014  seconds  and  approaches  a  steady 
value,  similar  to  the  final  steady  maximum  throttle  value  (not  the  quasi-steady  value), 
after  0.019  seconds  and  the  completion  of  the  analysis. 
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Thrust  vs.  Time 


Figure  6.38:  Transient  Thrust  Behavior  for  Inviscid  Wall  Dump  Combustor  Design 


Figure  6.39  shows  the  transient  behavior  of  the  total  entropy  at  the  exit  plane  of 
the  ramjet.  The  extreme  oscillatory  behavior  is  most  evident  in  this  case  with  the  dump 
combustor.  Similar  trends  are  observed  between  the  transient  thrust  and  entropy  values  at 
the  initiation  of  the  throttling  maneuvers  as  discussed  for  the  inviscid  wall  conventional 
design.  As  remarked  for  the  thrust  results  in  which  a  quasi-steady  thrust  value  was 
achieved  before  reaching  a  true  steady-state  value,  it  is  noted  here  that  the  entropy 


exhibits  the  same  behavior. 
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Entropy  vs.  Time 


Figure  6.39:  Transient  Total  Entropy  Behavior  for  Inviscid  Wall  Dump  Combustor 

Design 
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7.  AXISYMMETRIC  VISCOUS  WALL  RESULTS 

This  section  of  the  thesis  presents  results  and  analysis  for  the  conventional  (no¬ 
dump  combustor)  design  with  viscous  walls,  i.e.  when  utilizing  a  no-slip  boundary 
condition.  Only  the  conventional  (no-dump  combustor)  design  was  analyzed  in  this 
investigation  using  a  no-slip  boundary  condition  (viscous  wall).  This  was  done  due  to 
computational  and  time  limitations. 

The  viscous  wall  analysis  in  terms  of  approach  was  similar  to  the  inviscid  wall 
analysis  except  for  two  main  differences.  Firstly,  in  the  case  of  the  viscous  wall  analysis, 
a  tare  condition  was  not  established  before  the  first  maximum  throttle  condition  was 
established,  i.e.  the  analysis  used  the  maximumi  throttle  setting  as  the  initial  condition  for 
all  cases.  This  was  done  in  an  attempt  to  obtain  the  data  for  this  analysis  in  a  timely 
manner.  In  the  absence  of  hysteresis,  it  is  believed  that  it  should  have  little  effect  on  the 
values  obtained  from  the  subsequent  maximum  and  minimum  throttle  values  obtained 
from  CFD.  Secondly,  for  the  viscous  wall  analysis,  there  are  two  minimum  throttle 
settings  found.  This  is  the  result  of  establishing  a  new  minimum  throttle  setting  by 
increasing  heat  addition  due  to  the  absence  of  a  normal  shock  in  the  viscous  wall  analysis 
for  the  previous  minimum  throttle  setting.  The  two  minimum  throttle  settings  are 
therefore  referred  to  as  minimumi  (corresponding  to  the  previous  minimum)  and 
minimum2,  with  minimum2  being  the  actual  minimum  throttle  setting  based  upon  the 
definition  of  minimum  throttle  (normal  shock  residing  at  maximum  area  downstream  of 
inlet  and  upstream  of  heat  addition).  Figure  7.1  summarizes  the  theoretical  and  CFD 
based  maximum/minimum  throttle  settings  for  the  viscous  wall  analysis. 
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Figure  7.1:  Summary  of  Theoretical  and  CFD  Based  Maximum/Minimum  Throttle  Heat 

Additions  for  Viscous  Wall  Analysis 


7.1.  BLEED  BOUNDARY  CONDITION 

For  the  viscous  wall  analysis,  a  bleed  condition  was  specified  in  VULCAN  in  an 
attempt  to  prevent  the  ramjet  from  unstarting.  In  VULCAN,  the  bleed  boundary  condition 
was  set  by  specifying  the  amount  of  mass  flow  rate  to  be  extracted  from  the  mainstream 
flow  by  means  of  varying  the  back  pressure  across  the  bleed  region.  The  no-slip 
boundary  condition  imposed  created  a  boundary  layer  that  would  (in  the  presence  of  the 
adverse  pressure  gradient  in  the  ramjet)  separate  and  develop  into  a  recirculating  region 
just  upstream  of  heat  addition.  This  recirculation  region  would  then  grow  and  propagate 
upstream.  The  main  flow  would  see  this  recirculation  region  as  a  flow  blockage.  This 
would  compound  the  problem  and  would  lead  to  the  eventual  unstart  of  the  ramjet,  even 
at  minimum  throttle.  In  order  to  control  this  problem,  bleed  was  attempted.  Initially,  5% 
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of  the  air  mass  flow  rate  was  specified  to  be  bled,  (corresponding  to  5%  bleed)  i.e. 
removed  at  the  bleed  boundary.  It  was  determined  that  this  was  insufficient  bleed,  i.e.  the 
ramjet  would  still  unstart.  It  was  found  that  10%  bleed  proved  to  be  effective  in 
preventing  the  ramjet  from  unstarting;  however,  the  bleed  had  obvious  detrimental  effects 
on  the  thrust  performance  due  to  the  dependence  of  thrust  on  air  mass  flow  rate  at  the  exit 
of  the  engine.  Figure  7.2  shows  the  location  of  the  ramjet  wall  where  the  bleed  condition 
was  applied. 


Bleed  Boundary 
Condition 
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■  " 
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Axial  Location  [m] 


Figure  7.2:  Location  of  Bleed  Boundary  Condition  for  Viscous  Wall  Analysis 


Figure  7.3  shows  a  detail  of  the  bleed  region  including  the  stream  lines.  It  can  be 
seen  that  the  normal  shock  is  arrested  just  aft  of  the  bleed  boundary.  An  increase  in  the 
Mach  number  is  observed  just  before  the  normal  shock  at  the  bleed  region  due  to  the  flow 
expanding  in  order  to  compensate  for  the  loss  of  mass  flow  rate  out  of  the  ramjet. 
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Figure  7.3:  Mach  Number  Contours  and  Stream  Lines  of  Zoom  in  Portion  Where  Bleed 
Boundary  Condition  is  Specified  for  Viscous  Wall  Analysis 


7.2.  STEADY-STATE  THROTTLING 

Figure  7.4  shows  the  L2  of  the  residual  as  well  as  the  locations  of  heat  change  for 
the  viscous  wall  conventional  design  (no-dump).  This  information  is  also  summarized  in 


Table  7.1. 
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Figure  7.4:  L2  of  the  Residual  vs.  Iteration  for  Viscous  Wall  Conventional  Design 


Table  7.1:  Summary  of  Heat  Addition  and  Iteration  Count  for  Viscous  Conventional 

Design  Analysis 


Maximum 

Throttle 

Minimum  1 
Throttle 

Minimum! 

Throttle 

Maximum 

Throttle 

Heat  Addition  [kW] 

180.00 

0.72 

54.00 

180.00 

Current  Iterations 

500,000 

500,000 

700,000 

400,000 

Total  Iteration  Count 

500,000 

1,000,000 

1,700,000 

2,100,000 
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7.2.1.  Fluid  Dynamics.  Figure  7.5  shows  the  steady-state  Mach  number 
contours  achieved  for  the  viscous  wall  analysis  of  the  conventional  design  with  10% 
bleed.  The  normal  shock  is  just  aft  of  the  bleed  region  for  both  the  maximum  throttle 
cases.  The  normal  shock  has  degraded  to  an  oblique  shock  at  the  end  of  the  constant  area 
portion  of  the  ramjet  for  the  minimumi  throttle  setting.  Heat  rate  was  then  increased, 
causing  the  oblique  shock  from  the  minimumi  throttle  case  to  propagate  upstream.  This 
oblique  shock  propagates  upstream  but  still  remains  downstream  of  the  bleed  region  and 
turns  into  a  normal  shock  near  the  center-line  of  the  ramjet  for  the  minimumi  throttle 
setting.  The  shock  system  locations  for  the  viscous  wall  analysis  were  found  to  be  very 
sensitive  to  heat  addition  and  therefore  a  conventional  minimum  throttle  where  a  normal 
shock  was  established  within  the  constant  area  portion  of  the  ramjet  was  unattainable. 
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Figure  7.5:  Viscous  Wall  Summary  of  Mach  Number  Contours  at  Major  Throttle  Settings 

for  Conventional  Design 
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7.2.2.  Performance  and  Entropy  Results  and  Analysis.  The  thrust  data 
achieved  for  the  steady-state  analysis  of  the  conventional  ramjet  with  viscous  walls  and 
bleed  is  contained  in  Table  7.2. 


Table  7.2:  Summary  of  Thrust  Values  for  Steady-State  Viscous  Wall  Conventional 

Design 


Maximumi 

Throttle 

Minimumi 

Throttle 

Minimumi 

Throttle 

Maximumi 

Throttle 

Thrust  [N] 

33.6 

-23.8 

-15.9 

33.9 

Figures  7.6,  7.7,  7.8  and  7.9  provide  breakdowns  of  the  entropy  per  mass  for  the 
viscous  wall  ramjet  results  in  terms  of  the  various  loss  mechanisms  for  the  maximumi, 
minimumi,  minimurm,  and  maximumi  throttle  settings  respectively.  The  entropy  shown 
after  the  point  where  the  bleed  condition  has  been  specified  has  been  corrected  due  to  the 
bleed  boundary  condition  in  order  to  maintain  a  positive  value  of  entropy  at  the  exit  of 
the  ramjet,  and  also  for  comparative  purposes  (bleed  removes  entropy  associated  with  the 
removed  mass).  This  was  done  by  using  as  a  reference  the  entropy  per  mass  just  after  the 
bleed  region,  instead  of  at  the  inlet  face,  and  then  redefining  the  total  entropy  generated  at 
that  point  (just  aft  of  bleed)  to  be  equal  to  the  entropy  generated  from  friction  and 
convective  heat  transfer  in  the  inlet  (up  to  the  location  just  after  bleed).  This  means  that 
the  entropy  at  the  inlet  face  is  no  longer  zero,  since  the  total  entropy  is  reset  at  the  bleed 
location  to  account  for  the  removed  mass  and  entropy,  but  accounts  for  upstream  friction 
and  convective  heat  transfer  in  the  inlet.  This  however  assumes  that  there  is  no  entropy 
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generated  due  to  shocks  from  the  inlet  to  the  end  of  the  bleed  region.  This  is  a  relatively 
safe  assumption  after  inspection  of  Figures  7.6  through  7.9  (entropy  due  to  shocks  in  that 
location  is  nearly  zero). 
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Figure  7.6:  Entropy  Breakdown  for  Viscous  Wall  Conventional  Design  Maximumi 
Throttle  Case  (note  shifted  entropy  references  due  to  bleed) 


Figure  7.6  shows  the  effect  of  referencing  the  entropy  throughout  the  ramjet  to  the 
entropy  just  aft  of  the  bleed  region.  The  amount  of  entropy  lost  associated  with  the  mass 
that  is  extracted  across  the  bleed  boundary  condition  is  comparable  to  the  amount  of 
entropy  generated  downstream;  this  is  true  for  all  throttling  cases  presented  in  the  viscous 
wall  analysis.  Since  10%  of  the  mass  flow  rate  is  being  extracted,  approximately  a  10% 
decrease  in  the  total  entropy  entering  the  ramjet  is  observed  (slightly  more  than  10%  of 
the  entering  entropy  is  extracted  due  to  the  higher  concentration  of  entropy  within  the 
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boundary  layer  that  is  extracted  across  the  bleed  boundary).  Also  note  that  the  amount  of 
entropy  generated  throughout  the  ramjet  is  similar  to  the  amount  lost  across  the  bleed 
boundary  condition. 
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Figure  7.7:  Entropy  Breakdown  for  Viscous  Wall  Conventional  Design  Minimumi 

Throttle  Case 


Figure  7.7  shows  the  entropy  breakdown  for  the  minimumi  throttle  case.  This  plot 
exhibits  unsteadiness,  evidenced  by  the  decrease  in  total  entropy  and  thus  entropy 
generation  due  to  shocks  across  the  exit  nozzle.  For  this  minimumi  throttle  case,  more 
entropy  is  bled  from  the  system  than  is  generated  throughout  the  ramjet.  As  would  be 
expected  for  the  viscous  wall  (with  large  gradients  produces  by  the  no-slip  boundary 
condition),  entropy  generation  due  to  friction  is  much  larger  then  seen  for  the  inviscid 
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wall  analysis.  Note  the  absence  of  any  spike  in  entropy  due  to  the  absence  of  a  normal 
shock. 
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Figure  7.8:  Entropy  Breakdown  for  Viscous  Wall  Conventional  Design  Minimurm 

Throttle  Case 


The  presence  of  a  pronounced  shock  system  (oblique  shock  turning  into  a  normal 
shock  near  the  center- line)  is  evident  in  Figure  7.8  (minimum 2  throttle  condition)  at 
approximately  0.11  meters  from  the  ramjet  inlet.  It  is  interesting  to  note  that  the  shock 
system  is  not  actually  in  direct  contact  with  the  region  of  the  bleed  condition.  For  this 
new  minimum  throttle  case,  although  the  heat  added  (54  kW)  is  75  times  larger  than  in 
the  inviscid  wall  conventional  design  minimum  throttle  case  (0.72  kW)  the  entropy  due  to 
friction  is  similar  to  the  entropy  due  to  ideal  heat  transfer. 
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Total  Entropy  Flow  Rate  vs.  Axial  Location 
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Figure  7.9:  Entropy  Breakdown  for  Viscous  Wall  Conventional  Design  Maximuni2 

Throttle  Case 


Figure  7.9  provides  the  entropy  breakdown  for  the  viscous  wall  case  (maximum2 
throttle  condition).  It  shows  similar  characteristics  as  discussed  for  Figure  7.6. 

Figures  7.10,  7.11,  7.12,  and  7.13  show  the  percentage  breakdown  of  entropy  at 
the  exit  plane  of  the  ramjet  for  the  maximumi,  minimumi,  minimurm,  and  maximuni2 


throttle  settings,  respectively. 


86 


Figure  7.10:  Percentage  Breakdown  of  Entropy  Due  to  its  Mechanisms  for  Maximumi 
Throttle  Case  at  the  Exit  Plane  on  Viscous  Wall  Conventional  Design 


For  the  maximumi  throttle  case  (Figure  7.10),  it  can  be  seen  that  68.4%  of  the 
entropy  generation  is  due  to  ideal  heat  transfer.  The  combined  effect  of  entropy 
generation  due  to  friction  and  convective  heat  transfer  now  accounts  for  18.8%  of  the 
total,  whereas  in  the  inviscid  wall  analysis  for  the  conventional  design  the  combined 
effect  of  entropy  generation  due  to  friction  and  convective  heat  transfer  only  accounted 


for  1.5%  of  the  total  entropy. 
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Figure  7.11:  Percentage  Breakdown  of  Entropy  Due  to  its  Mechanisms  for  Minimumi 
Throttle  Case  at  the  Exit  Plane  on  Inviscid  Wall  Conventional  Design 


Figure  7.11  (entropy  results  at  exit  plane  for  minimumi  throttle  case)  shows  that 
while  the  impact  of  entropy  generation  due  to  shocks  decreases  from  76.8%  in  the 
inviscid  wall  conventional  design  analysis  to  56.5%  (as  previously  obtained)  for  the 
viscous  wall  conventional  design  analysis,  the  entropy  due  to  friction  has  increased  from 


0.2%  to  28.8%. 
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Viscous  Wall  Entropy  Allocation,  Minimum2 

13.7% 


Figure  7.12:  Percentage  Breakdown  of  Entropy  Due  to  its  Mechanisms  for  Minimumi 
Throttle  Case  at  the  Exit  Plane  on  Viscous  Wall  Conventional  Design 


Figure  7.12  shows  the  breakdown  of  entropy  generation  due  to  loss  mechanisms 
for  the  minimum 2  throttle  setting.  This  throttle  setting  increased  the  impact  of  entropy 
generated  from  ideal  and  non-ideal  heat  transfer  from  2.1%  to  46.4%,  as  compared  to  the 


minimumi  throttle  setting. 
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Viscous  Wall  Entropy  Allocation,  Maximum2 


Figure  7.13:  Percentage  Breakdown  of  Entropy  Due  to  its  Mechanisms  for  Maximumi 
Throttle  Case  at  the  Exit  Plane  on  Viscous  Wall  Conventional  Design 


The  maximumi  throttle  case  (Figure  7.13),  as  expected,  is  very  similar  to  the 
maximumi  throttle  case,  (Figure  7.10). 

A  summary  of  the  entropy  generation  allocation  at  the  exit  plane  of  the  ramjet  for 
all  the  throttling  cases  can  be  found  in  Table  7.3  for  the  viscous  wall  conventional  design. 
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Table  7.3:  Summary  of  Entropy  Allocation  at  Exit  Plane  of  Ramjet  for  All  Steady-State 


Throttling  Cases  for  Viscous  Wall  Conventional  Design 


Mechanism  [J/kg-K] 

Maximumi 

Minimum  i 

Minimurm 

Maximurm 

Friction 

21.45 

37.25 

59.34 

21.16 

Convective  HT 

31.63 

16.30 

6.83 

32.48 

Non-ideal  HT 

5.59 

1.57 

4.52 

5.62 

Ideal  HT 

192.84 

1.15 

72.38 

193.73 

Shocks 

30.52 

73.00 

22.68 

16.15 

Total 

282.04 

129.27 

165.74 

269.14 

7.3.  TRANSIENT  THROTTLING 

This  section  describes  results  for  transient  throttling  of  the  viscous  wall 
configuration  of  the  conventional  (no-dump)  ramjet. 

7.3.1.  Fluid  Dynamics.  The  transients  for  the  step  down  maneuver  can  be 
seen  in  the  snap  shots  of  the  Mach  number  contours  captured  in  Figure  7.14.  Each  frame 
represents  a  time  step  of  0.0010  seconds  for  a  total  elapsed  time  of  0.0120  seconds, 
approximately  covering  the  duration  of  the  transient  maneuver  and  achieving  the  new 
minimum  throttle  setting.  It  is  important  to  note  the  minimum  throttle  change  at  0.0050 


seconds. 


91 


■ 

Mach  Number:  0  0.5  1  1.5  ^ 

2.5  3 

t  =  0.0000  s 
t  =  0.0010  s 
t  =  0.0020  s 
t  =  0.0030  s 
t  =  0.0040  s 
t  =  0.0050  s 
t  =  0.0060  s 
t  =  0.0070  s 
t  =  0.0080  s 
t  =  0.0090  s 
t  =  0.0100  s 
t  =  0.0110  s 
t  =  0.0120  s 

0  0.05  0.1  0.15  0.2  0.25 

Axial  Location  [m] 

Figure  7.14:  Mach  Number  Contours  of  Step  Down  Maneuver  for  Viscous  Wall 

Conventional  Design 


Figure  7.15  shows  the  mass  flow  rate  through  the  ramjet  corresponding  to  Figure 
7.14.  The  location  of  the  bleed  boundary  location  is  apparent  by  observing  the  10% 
decrease  in  mass  flow  rate  across  this  region  at  the  inlet  throat.  The  no-slip  boundary 
condition  adds  a  considerable  amount  of  inherent  unsteadiness  to  the  system  as  can  be 
seen  by  the  transient  behavior  of  the  mass  flow  rate  in  Figure  7.15. 
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m  =  0.284  kg/s 


Figure  7.15:  Mass  Flow  Rate  Through  Ramjet  During  Transient  Step  Down  Maneuver 

for  Viscous  Wall  Conventional  Design 


The  transients  for  the  step  up  maneuver  can  be  seen  in  the  snap  shots  of  the  Mach 
number  contours  captured  in  Figure  7.16.  Each  frame  represents  a  time  step  of  0.0010 
seconds  for  a  total  elapsed  time  of  0.0040  seconds,  approximately  covering  the  duration 
of  the  transient  maneuver.  A  steady-state  solution  seems  to  be  reached  soon  relative  to 
the  step  down  maneuver,  which  was  similar  to  what  was  observed  for  the  inviscid  wall 
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conventional  design  analysis.  One  reason  for  this  could  be  the  relative  similarity  of  the 
shock  system  for  the  minimum2  and  maximum^  throttle  settings. 


M 

Mach  Number:  0  0.5 

1.5  ; 

2.5  3 

t  =  0.0000  s 
t  =  0.0010  s 
t  =  0.0020  s 
t  =  0.0030  s 
t  =  0.0040  s 

0  0.05  0.1  0.15  0.2  0.25 

Axial  Location  [m] 

Figure  7.16:  Mach  Number  Contours  of  Step  Up  Maneuver  for  Viscous  Wall 

Conventional  Design 


Figure  7.17  shows  the  mass  flow  rate  through  the  ramjet  corresponding  to  Figure 
7.16.  The  most  extreme  oscillatory  behavior  happens  directly  after  the  initiation  of  the 
step  up  maneuver  and  quickly  settles  to  a  relatively  steady-state  value.  The  decrease  in 
mass  flow  downstream  of  the  normal  shock  propagating  up  stream  is  only  slightly 
observable  since  the  shock  system  travels  upstream  only  a  relatively  short  distance. 
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m  =  0.284  kg/s 


Figure  7.17:  Mass  Flow  Rate  Through  Ramjet  During  Transient  Step  Up  Maneuver  for 

Viscous  Wall  Conventional  Design 
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7.3.2.  Performance  and  Entropy  Results  and  Analysis.  Figure  7.18  shows  the 
transient  thrust  behavior  for  the  viscous  wall  conventional  design  analysis.  The  trends 
observed  are  similar  in  behavior  to  the  inviscid  wall  conventional  design  analysis  and 
will  therefore  not  be  detailed.  The  major  differences  are  the  increased  oscillatory 
behavior  attributed  to  the  no-slip  boundary  condition  and  the  overall  decrease  in  the 
magnitude  of  the  thrust  values. 


Thrust  vs.  Time 


Figure  7.18:  Transient  Thrust  Behavior  for  Viscous  Wall  Conventional  Design 


Figure  7.19  shows  the  transient  behavior  of  the  total  entropy  at  the  exit  plane  of 
the  ramjet.  The  trends  observed  are  also  similar  in  behavior  to  the  inviscid  wall 
conventional  design  analysis  and  will  therefore  not  be  detailed.  The  major  difference  is 
the  increased  oscillatory  behavior  attributed  to  the  no-slip  boundary  condition.  It  is 
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interesting  to  note  that  the  magnitude  of  the  steady-state  total  entropy  at  the  exit  plane  of 
the  ramjet  for  the  viscous  wall  conventional  design  analysis  is  comparable  to  that  of  the 
inviscid  wall  conventional  design  analysis. 


Entropy  vs.  Time 


Time  [s] 

Figure  7.19:  Transient  Total  Entropy  Behavior  for  Viscous  Wall  Conventional  Design 
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8.  SUMMARY  AND  CONCLUSIONS 

VULCAN  CFD  code  was  used  to  analyze  the  steady-state  fluid  dynamics  and 
transient  entropy/performance  characteristics  of  ramjets  operating  at  a  flight  Mach 
number  of  3.0  and  standard  altitude  of  10  km.  Three  ramjets  were  analyzed,  an  inviscid 
wall  conventional  design,  an  inviscid  wall  dump  combustor  design  (to  prevent  unstart), 
and  a  viscous  wall  conventional  design  (utilizing  a  bleed  boundary  condition  to  prevent 
unstart).  The  formulation  of  the  overall  geometry  and  grids  of  the  ramjets  was  discussed. 
The  details  of  a  quasi-one-dimensional  analytical  model  used  to  determine  the  maximum 
and  minimum  throttle  settings  were  presented.  A  grid  convergence  study  was  then 
conducted  and  overall  grid  spacing  was  chosen  and  used  for  the  entirety  of  the 
investigation.  A  case  study  was  presented  to  determine  the  impact  of  two  different 
turbulence  models  on  the  flow-field  physics  with  focus  on  unstart  behavior.  Little 
variation  between  two  different  turbulence  models  in  terms  of  unstarting  behavior  was 
observed.  An  inviscid  wall  analysis  of  the  conventional  ramjet  was  then  conducted.  First, 
the  conventional  design  was  analyzed  as  a  basis  for  comparison.  Next,  the  dump 
combustor  was  analyzed.  Lastly,  a  viscous  wall  analysis  was  then  presented  utilizing  the 
conventional  design  only. 

The  inviscid  wall  analysis  demonstrated  that  the  dump  combustor  design  proved 
to  be  greatly  effective  at  preventing  unstart.  However,  the  dump  combustor  design  tended 
to  increase  the  entropy  generated  due  to  friction,  convective  heat  transfer,  and  non-ideal 
heat  transfer,  thereby  increasing  the  total  entropy  generated  within  the  ramjet  for  the 
maximum  throttle  setting.  At  minimum  throttle,  the  dump  combustor  had  a  lower  value 
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of  total  entropy  when  compared  to  the  conventional  design;  however,  this  is  most  likely 
due  to  the  absence  of  a  normal  shock  within  the  ramjet.  The  nature  of  the  dump 
combustor  also  added  to  inherent  unsteadiness  of  the  system. 

The  viscous  analysis  on  the  conventional  (no-dump)  configuration  utilized  a  bleed 
boundary  condition  at  the  inlet  throat  where  mass  flow  would  be  extracted  from  the  main 
flow.  10%  bleed  at  the  inlet  throat  was  found  effective  to  avert  unstart.  The  entropy  that 
was  bled  off  (associated  with  the  mass  flow  bled)  was  found  to  be  slightly  more  than  the 
total  entropy  entering  the  ramjet.  This  was  due  to  the  concentration  of  higher  entropy 
near  the  wall  upstream  of  the  inlet  throat  due  to  the  no-slip  boundary  condition  which 
was  then  subsequently  extracted  at  the  bleed  region.  It  was  observed  that  the  entropy 
generated  due  to  friction  and  convective  heat  transfer  were  the  only  mechanisms  of 
entropy  generation  that  significantly  increased  (the  dump  combustor  also  increase  non¬ 
ideal  heat  transfer).  The  no-slip  boundary  condition,  while  it  is  a  better  model  of  the  flow 
physics,  also  added  inherent  unsteadiness  to  the  system. 

The  transient  analysis  of  all  three  cases  studied  represented  the  most  extreme  (yet 
unphysical)  scenario  where  heat  was  added  instantaneously  between  throttle  maneuvers. 
The  transients  observed  were  relatively  short  lived  considering  the  time  frame  of  actual 
throttle  maneuvers  and  exhibited  vast  oscillatory  behavior.  If  the  assumption  of  steady- 
state  flow-field  variables  were  made  for  the  analysis  and  optimization  of  ramjet 
performance,  the  results  could  be  skewed  by  the  accumulation  of  these  deviations  from 
steady-state  values  over  time. 

While  there  was  no  clear  preferred  method  that  prevented  unstart  while  having  a 
minimal  impact  on  entropy  generation,  it  is  important  to  note  that  the  dump  combustor 
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(along  with  the  slip  boundary  condition  making  it  an  inviscid  wall  analysis)  did  not 
model  the  flow  physics  as  well  as  the  viscous  wall  analysis  on  the  conventional  design. 
However,  both  methods  were  highly  effective  at  preventing  unstart  and  exhibited 
expected  entropy  generation  characteristics. 


APPENDIX  A. 

SAMPLE  VULCAN  INPUT  DECK 
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t~)'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

$****************  Begining  of  general  control  data  ************* 

■k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


$ - Parallel  processing  control  data - 

PROCESSORS  24.0  (No.  of  cpus  to  use) 

MESSAGE  MODE  0.0  (Message  passing  strategy:  0=stnd., 


$ 

$ 

$ 

$ 


l=buf f ered) 

$ - Geometric  model  type - $ 

AXISYM  1.0  (twod,  axisym,  threed) 

$ - Grid  file  data - $ 

GRID  FORMAT  3.0  (l=s.b.form,  2=s.b.bin.,  3=m.b.form., 


4=m. b . bin . ) 

GRID  0.0  ( 0=plot3d->3d  ;  plot2d->2d/axi , 

l=plot3d->all) 

double  grid  vl9  AWALL.grd 

GRID  SCALING  FACTOR  0.25  (Converts  grid  units  to  meters) 

$ - Restart  file  data - $ 


RESTART  IN  1.0 

vl9  AWALL  Q2500  2  4. restart 
RESTART  OUT  1.0 

vl9  AWALL  Q2500  2_5 . restart 
RESTART  OUT  INTERVAL  10000.0 


$ -  Output 

WARNING  MESSAGES  0.0 

3=both  ) 

PLOT  ON  3.0 

4=m. b . unf rm. ) 

PLOT  NODES  0 . 0 

averaged  to  the  nodes) 
PLOT  FUNCTION  9.0 


control  data - $ 

(0=none,  l=wall  funct.,  2=temp.  limit, 

(l=s.b.frm.,  2=s . b . unf rm. ,  3=m.b.frm., 

(Create  PLOT3D  files  using  data 

(Create  PLOT3D  function  file 


containing  variables  below) 


DENSITY 
VELOCITY 
PRESSURE 
TEMPERATURE 
MACH  NO. 

LAM.  VIS. 

EDDY  VIS.  RATIO 
TOTAL  ENTH. 

GAMMA 

OUTPUT  TIME  HISTORY  1000.0  (iteration  interval  between  time 
history  writes) 

timhis  vl9  AWALL  Q2500  2  5.tim 

$ -  Gas  thermo,  diffusion,  and  reaction  model  data  --$ 

GAS/THERMO  MODEL  1.0  (0=CPG,  1=TPG,  2=n/a) 

CHEMISTRY  MODEL  0.0  (0=frozen,  l=finite  rate,  2=n/a) 

GLOBAL  VISCOUS  0.0  (solve  the  Navier-Stokes  equations) 

$ - Transport  model  data - $ 

VISCOSITY  MODEL  1.0  ( l=Sutherlands  law) 

CONDUCTIVITY  MODEL  0.0  (0=Prandtl  no.,  l=Wassilej's  law) 

UNIV.  GAS  CONST.  8314.34 
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NO.  OF  CHEMICAL  SPECIES  2.0 

/ share /apps/Vul can /Ver_6 . 0 . 2 /Data_base/gas_mod . Lewis_3 
N2  02 

0.7686  0.2314 

$ - Reference  condtion  data - $ 

ANGLE  REF.  FRAME  0.0  (0=alpha  in  xy  plane,  l=alpha  in  xz 

plane) 

ALPHA  0.0  (angle  of  attack  measured  C.C.W  in 

degrees ) 

NONDIM  1.0  ( 0=non . dimen . ,  l=dimen.  static, 

2=dimen.  total) 

MACH  NO .  3.0 

STATIC  TEMP.  223.15 

STATIC  PRESS.  26436.267594 

LAM.  PRANDTL  NO.  0.72 

LAM.  SCHMIDT  NO.  0.22 

TURB.  PRANDTL  NO.  0.90 

TURB.  SCHMIDT  NO.  0.90 

$ -  Turbulence  Model  Data  - $ 

TURB.  MODEL 

MENTER-SST  (SPALART,  MENTER,  MENTER-SST,  K-OMEGA) 

TURB.  INTENSITY  0.05 

TURB.  VISC.  RATIO  0.10 

BOUSSINESQ  REY.  STRESS  0.0 

NO  2/3  RHOK  IN  REY.  STRESS  0.0 


$ -  Runge-Kutta  scheme  coefficients  - $ 

NSTAGE  3.0  (no.  of  Runge-Kutta  Stages) 

0.333333333333,  0.5,  1.0 

$ - Boundary  and  cut  control - $ 

FLOWBCS  50.0  (no.  of  boundary  conditions  to  be 

specified) 

BCGROUPS  28.0  (no.  of  boundary  condition  groups) 

CUTBCS  23.0  (no.  of  C(0)  conectivity  conditions 

to  be  specified) 

PATCHBCS  0.0  (no.  of  non-C(O)  conectivity 

conditions  to  be  specified) 

IGNITION  SUB-BLOCKS  2.0  (no.  of  ignition  sub-blocks) 

TIME  HISTORY  I/O  24.0  (no.  of  time  history  sub-blocks) 

BLOCKS  24.0  (no.  of  blocks) 

BLOCK  CONFIG.  24.0  (no.  of  lines  of  block  configurations 

input) 

BLK  I-VISCS  J-VISCS  K-VISCS  TURB  PLOT  SOLVER  REGION 
IT  T  N  Y  Y  E/A  1 

2  T  T  N  Y  Y  E/A  1 

3  T  T  N  Y  Y  E/A  1 

4  T  T  N  Y  Y  E/A  1 

5  T  T  N  Y  Y  E/A  1 

6  T  T  N  Y  Y  E/A  1 

7  T  T  N  Y  Y  E/A  1 

8  T  T  N  Y  Y  E/A  1 

T  T  N  Y  Y  E/A  1 


9 
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10  T  T  N  Y  Y  E/A  1 

11  T  T  N  Y  Y  E/A  1 

12  T  T  N  Y  Y  E/A  1 

13  T  T  N  Y  Y  E/A  1 

14  T  T  N  Y  Y  E/A  1 

15  T  T  N  Y  Y  E/A  1 

16  T  T  N  Y  Y  E/A  1 

17  T  T  N  Y  Y  E/A  1 

18  T  T  N  Y  Y  E/A  1 

19  T  T  N  Y  Y  E/A  1 

20  T  T  N  Y  Y  E/A  1 

21  T  T  N  Y  Y  E/A  1 

22  T  T  N  Y  Y  E/A  1 

23  T  T  N  Y  Y  E/A  1 

24  T  T  N  Y  Y  E/A  1 

REGION  CONFIG.  1.0  (no.  of  regions  the  blocks  are  grouped  into) 
$***************  Region  1  control  input  **********************$ 
ROE  KAPPA  LIMITER  LIM-COEF  ENTRP (U) 

ENTRP (U+a) 

3,  3,  3,  4,  4,  4,  2.0,  2.0,  2.0,  1.0,  1.0,  1.0,  1.0, 

1.0,  1.0 

FMGLVLS  NITSCG1  NITSCG2  NITSFG  # 1ST-ORD . -C . G . / ITER . 

RES. ; REL . , ABS. 

1  100000  0  -10.0  -10.0 
MG-CYCLE  COARSE  GRIDS  DQ-SMOOTH  DQ-CORR  DAMP-MEAN  DAMP-TURB 
I  0  0.25  0.50  1.0  0.5 

TURB  CONVECTION  DT  RATIO  NON-EQUIL  POINT-IMP  COMP  MODEL  CG 
WALL  BC 

1ST  0.1  25.0  N  N 

WMF 

SCHEME  TIME  STEP  IT-STATS  CFL-MIN  ADP-CFL  #CFL-VAL  VISC-DT 
IMP-BC  REG-REST 

DAF  SUBIT1  10  0.1  Y  2  Y  NY 

1  100000 

1.0  1.0 

TIME  STEP  SUB-ITS  RES-RED  TIME-ORDER  C-N  RELAX 
1.0e-8  5  -2.0  2ND  0.50 

i ***************  End  Qf  general  control  data  *******************; 

BC  GROUPS:  NAME  TYPE  OPTION 

inflow  REFFIX  PHYSICAL 

SUBMDOTO  SUBMDOTO  PHYSICAL 

Mass  Flow  Rate  Relax  Back  Pres 
0.0284404  0.5  -1.0 

top _ 1  AWALLM  PHYSICAL 

top_2  AWALLM  PHYSICAL 

top_3  AWALLM  PHYSICAL 

top_4  AWALLM  PHYSICAL 

top_5  AWALLM  PHYSICAL 

top_6  AWALLM  PHYSICAL 

top_7  AWALLM  PHYSICAL 
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top_8 

top_9 

top_10 

top_ll 

top_12 

top_13 

top_14 

top_15 

top_l 6 

top_17 

top_18 

top_l 9 

top_20 

top_21 

top  22 

top_23 

top_24 

bottom 

outflow 


AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AWALLM 

AXICL 

EXTRAP2 


PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 

PHYSICAL 


BC  NAME 

BLK 

FACE 

PLACE 

DIREC1 

BEGIN 

END  DIREC2 

BEGIN  END  IN- 

•ORDER 

inflow 

1 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

1 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  1 

1 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

2 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  2 

2 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

3 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  3 

3 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

4 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  4 

4 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

5 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

SUBMDOTC 

i  5 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

6 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  6 

6 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

7 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  7 

7 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

8 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  8 

8 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

9 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  9 

9 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

10 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  10 

10 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

11 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  11 

11 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

12 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  12 

12 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

13 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  13 

13 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

14 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  14 

14 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

15 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 
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top  15 

15 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

16 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  1 6 

16 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

17 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  17 

17 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

18 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  18 

18 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

19 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  1 9 

19 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

20 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  20 

20 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

21 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  21 

21 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

22 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  22 

22 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

23 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  23 

23 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

bottom 

24 

J 

MIN 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

top  24 

24 

J 

MAX 

I 

MIN 

MAX 

K 

MIN 

MAX 

0 

outflow 

24 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  NAME 

:  BLK 

FACE 

PLACE 

DIREC1 

BEGIN 

END 

DIREC2 

BEGIN 

END  INORDER 

CUT  1 

1 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  1 

2 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  2 

2 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  2 

3 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  3 

3 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT_3 

4 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  4 

4 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  4 

5 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  5 

5 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  5 

6 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  6 

6 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT_6 

7 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  7 

7 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  7 

8 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  8 

8 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  8 

9 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  9 

9 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  9 

10 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  10 

10 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  10 

11 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  11 

11 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  11 

12 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  12 

12 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  12 

13 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  13 

13 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT_13 

14 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  14 

14 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  14 

15 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  15 

15 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 
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CUT  15 

16 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  16 

16 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT_1 6 

17 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  17 

17 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  17 

18 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  18 

18 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  18 

19 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  19 

19 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  19 

20 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  2  0 

20 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  2  0 

21 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  21 

21 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  21 

22 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  22 

22 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  22 

23 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  23 

23 

I 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

CUT  23 

24 

I 

MIN 

J 

MIN 

MAX 

K 

MIN 

MAX 

0 

IGN  S.E 

i.  BLK  DIR1 

BEG  END 

DIR2 

BEG  END  DIR3  BEG 

END  T/P- 

-IGNITE 

SPARK1 

19 

I 

MIN  MAX 

J 

MIN  MAX  K 

MIN 

MAX 

-1250 

SPARK2 

20 

I 

MIN  MAX 

J 

MIN  MAX  K 

MIN 

MAX 

-1250 

TIM  HIS 

BLK 

DIRl  BEG 

END 

DIR2 

BEG 

END 

DIR3 

BEG 

END 

TIMHIS 

1  1 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

2  2 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

3  3 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

4  4 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

5  5 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

6  6 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

7  7 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

8  8 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

9  9 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

10  10 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

11  11 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

12  12 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

13  13 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

14  14 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

15  15 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

16  16 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

17  17 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

18  18 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

19  19 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

20  20 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

21  21 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

22  22 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

23  23 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

TIMHIS 

2  4  2  4 

I 

MIN 

MAX 

J 

MIN 

MAX 

K 

MIN 

MAX 

APPENDIX  B. 

POST  PROCESSING  MATLAB  CODE 
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function  timhis_plot3d  animate 
tic 

clear  all 
close  all 
clc 


S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-S-1 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOi' 

%  Author:  Jonathan  Sheldon 


This  program  takes  individual 
block  from  VULCAN  and  assembl 
using  those  solution  files,  s 
for  ease  of  later  calculation 
are  calc_parameters  and  corre 
parameters  and  the  name  file 
parameters . 


time  history  solution  files  for  each 
es  them  into  one  solution  file.  Then 
tores  each  parameter  in  its  own  matrix 
s  for  each  block.  Things  to  be  changed 
sponding  calculations  for  the  desired 
generated  to  properly  reflect  the 


n  blocks  =  24; 
n  timhis  =  100; 
n_parameters  =  10; 
calc_parameters  =  8; 

nn_parameters  =  n_parameters  +  calc_parameters ; 


parameter^ 

parametep 

parameter 

parameter^ 

parameter 

parameter^ 

parameter 

parameter 

parameter^ 

parameter 

parameter 


name 
name 
name 
name 
name 
name 
name 
name ( 8 , 
name ( 9 , 
name (10, 
name (11, 


(1, 

(2, 

(3, 

(4, 

(5, 

(6, 

(7, 


'Static  Density 
'u  Velocity 
'v  Velocity 
'Static  Pressure 
'Static  Temperature 
'Mach  Number 
'Laminar  Viscosity 
'Eddy  Vise.  Ratio 
'Total  Enthalpy 
' Gamma 


)  =  'Entropy 


parameter 

name ( 12 , 

)  =  'Entropy 

Gen . 

(Total) 

1  . 
r 

parameter 

name (13, 

)  =  'Entropy 

Gen . 

(friction) 

l  . 
r 

parameter 

name (14, 

)  =  'Entropy 

Gen . 

(conv.  HT) 

i  . 
r 

parameter 

name (15, 

)  =  'Entropy 

Gen . 

(non-ideal  HT) 

l  . 
r 

parameter 

name (16, 

)  =  'Entropy 

Gen . 

(ideal  HT) 

i  . 
r 

parameter 

name ( 17 , 

)  =  'Entropy 

Gen . 

(HT) 

l  . 
r 

parameter 

name (18, 

)  =  'Entropy 

Gen . 

( shocks ) 

l  . 
r 

Qdot  =  2500*72; 


r_i  =  0.0625/4; 
L  b  =  0.10/4; 
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% -  Create  timhis  File?  - 

% - Yes  =  1  - No  =  0 

Create  timhis  File  =  0; 


% -  Create  timhis  solution  File? 

% - yes  =  1  - No  =  0  --- 

Create  timhis  solution  File  =  0; 


% -  Create  Contour  Animation? 

% - Yes  =  1  - No  =  0 

Create  Contour  Animation  =  0; 


Create  CL  Plot  =  0; 


Create  Center-Line  Plot? 
Yes  =  1  - No  =  0 


Create  LG  movie  =  1; 


Create  Line  Graph  Movie/Entropy  Results? 
Yes  =  1  - No  =  0  - 


grid  =  dlmread /double  grid  vl9  AWALL.grd'); 

n  x  nodes  =  grid (2,1); 
n_y_nodes  =  grid (2,2); 
n  t  nodes  =  n  x  nodes*n  y_nodes; 


timhis  File 


if  Create  timhis_File  ==  1 
for  t  =  l:n  timhis 

write  label  =  [ ' timhis  ' , int2str (t) , ' . f ' ] ; 

fprintf (' Start  Creation  of  File:  %s  \n', write  label) 

read  label  1  = 


timhis  vl9 

AWALL 

Q2500 

2 

5  .  tim 

1 

, int2str (t) 

,  '  .f ' 

]  ; 

read 

label 

2  = 

timhis  vl9 

awall" 

Q2500 

2 

5  .  tim 

2 

_1  ' 

, int2str (t) 

,  '  .  f  ' 

]  ; 

read 

label 

3  = 

timhis  vl9 

awall" 

Q2500 

2 

5  .  tim 

_3 

1  ' 

, int2str (t) 

,  '  -  f  ' 

]  ; 

read 

label 

4  = 

timhis  vl9 

awall" 

Q2500 

2 

5  .  tim 

_4 

, int2str (t) 

,  '  -  f  ' 

]  ; 

read 

label 

5  = 

timhis  vl9 

AWALL 

Q2500 

JZ 

5  .  tim 

_5_ 

_1  ' 

, int2str (t) 

,  '  -  f  ' 

]  ; 

read 

label 

6  = 

timhis  vl9 

awall’ 

Q2500 

2 

5  .  tim 

6 

_1  ' 

, int2str (t) 

,  '  .  f  ' 

]  ; 

read 

label 

7  = 

timhis  vl9 

awall" 

Q2500 

2 

5  .  tim 

_7_ 

_1_' 

, int2str (t) 

,  '  -  f  ' 

]  ; 

read 

label 

8  = 

timhis  vl9 

awall" 

Q2500 

2 

5  .  tim 

_8 

_1  ' 

, int2str (t) 

,  '  -  f  ' 

]  ; 

read 

label 

9  = 

timhis  vl9 

awall" 

Q2500 

2 

5  .  tim 

9 

1  ' 

, int2str (t) 

,  '  -  f  ' 

] ; 

[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 
[ ' timhis 


read  label  10  = 


_vl9_ 

AWALL 

Q2500 

2 

5  .  tim 

_10_ 

1_' 

, int2str (t) 

,  '  -f 

read 

label 

11  = 

_vl9" 

awall" 

Q2500 

2 

5  .  tim 

_11_ 

1_' 

, int2str (t) 

,  '  -f 

read 

label 

12  = 

vl9" 

awall’ 

Q2500 

2 

5  .  tim 

12 

1  ' 

, int2str (t) 

,  '  -f 

read 

label 

13  = 

vl9" 

awall’ 

Q2500 

2 

5  .  tim 

_13_ 

1  ' 

, int2str (t) 

r  '  .f 

read 

label 

14  = 

_vl9" 

AWALL 

Q2500 

2 

5  .  tim 

_14^ 

1_' 

, int2str (t) 

,  '  -f 

read 

label 

15  = 

_vl9" 

awall" 

Q2500 

2 

5  .  tim 

_15_ 

1_' 

, int2str (t) 

,  '  -f 

read 

label 

16  = 

_vi  9] 

awall" 

Q2500 

2 

5  .  tim 

1 6_ 

1_' 

, int2str (t) 

,  '  -f 

read 

label 

17  = 

vl9" 

awall’ 

Q2500 

2 

5  .  tim 

_17^ 

1  ' 

, int2str (t) 

,  '  -f 

read 

label 

18  = 

vl9" 

awall’ 

Q2500 

2 

5  .  tim 

/L8_ 

i_' 

, int2str (t) 

,  '  -f 

read 

label 

19  = 

_vi  9~ 

AWALL 

Q2500 

2 

5  .  tim 

1 9__ 

1_' 

, int2str (t) 

,  '  -f 

read 

label 

20  = 

_vl9" 

awall" 

Q2500 

2 

5  .  tim 

2  0_ 

1_' 

, int2str (t) 

,  '  .f 

read 

label 

21  = 

_vl9" 

awall" 

Q2500 

2 

5  .  tim 

2i^ 

1_' 

, int2str (t) 

,  '  -f 

read 

label 

22  = 

vl  9" 

awall’ 

Q2500 

_2 

5  .  tim 

_2  2^ 

1  ' 

, int2str (t) 

,  '  -f 

read 

label 

23  = 

vl  9" 

awall" 

Q2500 

2 

5  .  tim 

23_ 

1  ' 

, int2str (t) 

,  '  -f 

read 

label 

24  = 

_vi  9] 

AWALL 

Q2500 

2 

5  .  tim 

_2  4_ 

l_ ' 

, int2str (t) 

,  '  -f 

timhis_l  = 
timhis_2  = 
timhis_3  = 
timhis_4  = 
timhis_5  = 
timhis__6  = 
timhis_7  = 
timhis_8  = 
timhis_9  = 
timhis_10 
timhis_ll 
timhis_12 
timhis_13 
timhis_14 
timhis_15 
timhis_l 6 
timhis_17 
timhis_18 
timhis_l 9 
timhis_20 
timhis_21 
timhis_22 
timhis__23 
timhis  24 


dlmread(read  label_ 
dlmread(read  label_ 
dlmread(read  label_ 
dlmread(read  label_ 
dlmread (read_label_ 
dlmread(read  label_ 
dlmread (read  label_ 
dlmread (read  label_ 
dlmread (read  label_ 
=  dlmread (read_label 
=  dlmread (read_label 
=  dlmread (read_label 
=  dlmread (read_label 
=  dlmread (read_label 
=  dlmread (read_label 
=  dlmread (read  label 
=  dlmread (read_label 
=  dlmread (read_label 
=  dlmread (read  label 
=  dlmread (read_label 
=  dlmread (read_label 
=  dlmread (read_label 
=  dlmread (read^label 
=  dlmread (read  label 


1) 

2) 

3) 

4) 

5) 

6) 
7) 


10) 

11) 

12) 

13) 

14) 

15) 

16) 

17) 

18) 

19) 

20) 
21) 
22) 

23) 

24) 


RPB 


=  ceil (n_x_nodes*n_y_nodes/3) ; 
dlmwrite (write  label, n  blocks, '  ') 

dlmwrite (write  label, timhis_l (2 , 1 : 2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  '  1 

dlmwrite  (write  label ,  timhis__2  (2 , 1 :  2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  '  ' 

dlmwrite (write  label , timhis_3 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter', 
dlmwrite (write  label , timhis_4 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append '  ,  'delimiter', 
dlmwrite (write  label , timhis_5 (2 , 1 : 2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  ' 
dlmwrite (write  label ,  timhis__6  (2 , 1 :  2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  ' 
dlmwrite  (write  label ,  timhis__7  (2 , 1 :  2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  1  ' 

dlmwrite (write  label , timhis_8 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append '  ,  'delimiter', 
dlmwrite (write  label , timhis_9 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter', 
dlmwrite (write  label , timhis_10 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  '  ' 

dlmwrite (write  label, timhis_ll (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  '  ' 

dlmwrite (write  label , timhis_12 (2 , 1 : 2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  1  ' 

dlmwrite (write  label , timhis_13 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter', 
dlmwrite (write  label , timhis_14 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter', 
dlmwrite (write  label , timhis_15 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  ' 
dlmwrite (write  label , timhis_l 6 (2 , 1 : 2 )  ,  ' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  '  ' 

dlmwrite (write  label , timhis_17 (2 , 1 : 2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  1  ' 

dlmwrite (write  label , timhis_l 8 (2 , 1 : 2 ) ,  ' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append '  ,  'delimiter', 
dlmwrite (write  label , timhis_l 9 (2 , 1 : 2 ) , ' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append '  ,  'delimiter', 
dlmwrite (write  label , timhis_20 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  '  ' 

dlmwrite (write  label , timhis_2 1 (2 , 1 : 2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  '  ' 

dlmwrite (write  label , timhis_22 (2 , 1 : 2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter',  1  ' 

dlmwrite (write  label , timhis_23 (2 , 1 : 2 ),' -append ' ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter', 
dlmwrite (write  label , timhis_24 (2 , 1 : 2 ),' -append '  ,  'delimiter', 
dlmwrite (write  label , n_parameters ,' -append ' ,  'delimiter', 
dlmwrite (write  label , timhis_l (RPBt2 : end, 1:3), ' -append ' , 
delimiter',  '  ’,  'precision',  '%.15e') 

dlmwrite (write  label , timhis_2 (RPB+2 : end, 1:3),  -append 1  , 
delimiter',  1  ',  'precision',  '%.15e') 
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dlmwrite (write_label , timhis_3 (RPB+2 : end,  1:3),  ' -append '  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_4 (RPB+2 : end, 1:3), ' -append ' , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_5 (RPB+2 : end, 1:3),  ' -append ' , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_6 (RPB+2 : end, 1:3), ' -append ' , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_7 (RPB+2 : end, 1:3), ' -append ' , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_8 (RPB+2 : end, 1:3), ' -append ' , 
'delimiter',  '  ',  'precision  ,  '%.15e') 

dlmwrite (write_label , timhis_9 (RPB+2 : end, 1:3), ' -append ' , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_10 (RPB+2 :end, 1:3), ' -append  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_l 1 (RPB+2 : end, 1:3), ' -append  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_12 (RPB+2 : end, 1:3), ' -append  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_13 (RPB+2 : end, 1:3),  ' -append '  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite  (write_label ,  timhis___14  (RPB+2  :  end,  1:3),  '  -  append  '  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_15 (RPB+2 :end, 1:3), ' -append  , 
'delimiter',  '  ,  'precision',  '%.15e') 

dlmwrite (write_label , timhis_l 6 (RPB+2 : end, 1:3), ' -append  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_17 (RPB+2 :end, 1:3), ' -append  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_l 8 (RPB+2 : end, 1:3),  ' -append  1 , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_l 9 (RPB+2 : end,  1:3),  ' -append ' , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_20 (RPB+2 :end, 1:3), ' -append  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_2 1 (RPB+2 : end,  1:3),  ' -append  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_22 (RPB+2 : end, 1:3), ' -append  , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite (write_label , timhis_2  3 (RPB+2 : end, 1:3),  ' -append  1 , 
'delimiter',  '  ',  'precision',  '%.15e') 

dlmwrite  (write_label ,  timhis__2  4  (RPB+2  :  end,  1:3),  '  -  append  ’  , 
'delimiter',  '  ',  'precision',  '%.15e') 

end 
toe 

fprintf (' Creation  of  timhis  X.f  Files  Completed  \n  \n') 


end 


o\°  o\°  o\° 
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Solution  File 


if  Create  timhis_solution_File  ==  1 
for  t  =  l:n  timhis 

read_label  =  [ ' timhis_' , int2str (t) , ' . f ' ] ; 

write  label  =  ['timhis  ' , int2str (t) , '  solution . f '] ; 

fprintf (' Start  Creation  of  File:  %s  \n ' , write_label ) 

solution  =  dlmread (read^label) ; 

dlmwrite (write  label, n  blocks) 

for  n  =  l:n  blocks 

dlmwrite (write^label ,  [n  x  nodes, n_y  nodes],  '- 
append' , 'delimiter' , '  ' ) 

dlmwrite (write_label , nn_parameters , '-append', 'delimiter1, '  ') 

end 

for  b  =  l:n  blocks 

clear  parameter  1  parameter  2  parameter  3  parameter  4 
parameter  5  parameter_6  parameter  7  parameter_8 
parameter  9  parameter  10 
if  b  ==  1 

first_row  =  2*n^blocks+2 ; 
last_row  = 

ceil (n_parameters*n  x  nodes*n  y  nodes/3+f irst  row-1); 

else 

first  row  =  last  row+1; 
last_row  = 

first  row+ceil (n_parameters*n  x  nodes*n_y  nodes/3-1); 

end 

solution  block  =  solution ( first_row : last_row, 1 : 3 ) ; 

% -  Build  Matricies  of  Parameters  - 

parameter  =  zeros (n  y^nodes,n  x  nodes,n_parameters) ; 
n  =  1; 
m  =  2  ; 

1  =  1; 

for  p  =  1 : n_parameters 

for  node  =  1 : n_parameters*n  t  nodes 

if  node  >  (p-l)*n  t  nodes  &&  node  <  p*n  t  nodes+1 
if  m  ==  4 

n  =  n+1; 
m  =  1  ; 

end 

if  node  ==  l*n  x  nodes+1 
1  =  1+1; 

end 

i  =  node  -  n  x  nodes* (1-1) ; 
j  =  1  -  (p-l)*n_y_nodes; 

parameter ( j , i , p)  =  solution_block (n, m) ; 
m  =  m  +  1  ; 

end 

end 

end 


Start  Calculations 
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Entropy 


[kg_N2/kg_air] 

[kg_02/kg_air] 


mF  N2  =  0.7686;  %mass  fraction  of  N2 

mF  02  =  0.2314;  %mass  fraction  of  02 

MW  N2  =  28.016;  %molecular  weight  of  N2  [kg  N2/kmol_N2] 

MW  02  =  32.000;  %molecular  weight  of  02  [kg  02/kmol  02] 

kmoles  N2  =  mF  N2/MW  N2;  %kmoles  of  N2 
kmoles  02  =  mF_02/MW  02;  %kmoles  of  02 
X  N2  =  kmoles  N2/ (kmoles  N2+kmoles  02) ; 


R  univ 
R_N2  = 
R  02  = 


=  8314.34; 
R_univ/MW_ 
R  univ/MW 


[ kmo 1 e  s ] 

[ kmoles ] 

%mole  fraction  of  N2 
[kmol  N2/kmol_air] 
%mole  fraction  of 
02  [kmol_02/kmol  air] 
%universal  gas  constant  [J/kmol-K] 

N2;  %N2  gas  constant  [J/kg  N2-K] 

02;  %02  gas  constant  [J/kg_02-K] 


X  02  =  kmoles  02/ (kmoles  N2+kmoles  02) 


%Coef f icients  for  polynomial  curve  fit  for  C_P(T) 
%C_P(T)  =  a  +  b*T  +  c*TA2  +  d*TA3  +  e*T&4 

a_N2  =  0 . 31459E+01*R_N2; 
b_N2  =  0 . 99154E-03*R_N2; 
c_N2  =  -0 ,22912E-06*R_N2; 
d_N2  =  0 . 12181E-10*R_N2; 
e_N2  =  0 . 11024E-14*R_N2; 

a_02  =  0 . 30809E+01*R_02; 
b_02  =  0 . 16962E-02*R_O2; 
c__02  =  -0 . 76334E-06*RJ02; 
d_02  =  0 . 17140E-09*R_02; 
e_02  =  -0.14116E-13*RJ02; 

s_ref_N2  =  191 . 61*1000/MW_N2;  %[J/kg_N2-K] 
s_ref_02  =  205 . 142*1000/MW_02; % [ J/kg_02-K] 


P_ref  =  100000.0;  %[Pa] 
T  ref  =  298.0;  % [K] 


P  =  parameter (:,:, 4 ) ;  %[Pa] 

T  =  parameter ( : , : , 5) ;  % [K] 

s_N2  =  s_ref_N2  +  a_N2*log (T/T_ref )  +  b_N2* (T-T_ref )  +  ... 
l/2*c_N2* (T. A2-T_refA2)  +  l/3*d_N2* (T . A3-T_ref A3)  +  ... 
l/4*e_N2* (T. A4-T_refA4)  -  R_N2*log (P/P_ref )  - 
R_N2*log (X_N2) ; 

s  02  =  s  ref_02  +  a  02*log (T/T^ref )  +  b_02* (T-T  ref)  +  ... 
l/2*c_02* (T. A2-T_refA2)  +  l/3*d_02*7T . A3-T_ref A3)  +  ... 
l/4*e_02* (T. A4-T_refA4)  -  R_02*log (P/P_ref )  - 
R_02*log (X_02 ) ; 

s  air  =  mF  N2*s  N2  +  mF  02 *s  02; 
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parameter 11 )  =  s_air; 

% -  Grid  Construction  Per  Block' 


if  b  ==  1 

first  row  x 
last_row_x 
first  row  y 
last_row_y 

else 

first  row  x 
last_row  x 
first  row^y 
last_row_y 

end 


n  blocks  +  2; 
first  row  x  + 
last  row  x  + 
first  row  y  + 

last  row_y  + 
first  row  x  + 
last  row  x  + 
first  row  y  + 


(n_y_nodes-l ) ; 
1; 

(n_y_nodes-l ) ; 
1; 

(n_y_nodes-l ) ; 
1; 

(n_y_nodes-l ) ; 


x  grid  =  0 . 25*flipud (grid (first  row  x:last  row  x,  :  )  )  ; 
y_grid  =  0 . 25*flipud (grid (f irst_row_y : last_row_y,  : ) )  ; 


Entropy  Break  Down 


rho 

=  flipud (parameter ( 

r 

,  l ) )  ; 

u 

=  flipud (parameter ( 

r 

,2)); 

V 

=  flipud (parameter ( 

r 

,3)); 

P 

=  flipud (parameter ( 

r 

,4)); 

T 

=  flipud (parameter ( 

r 

,5)); 

M 

=  flipud (parameter ( 

f 

,  6 )  )  ; 

mu  lam 

=  flipud (parameter ( 

t 

,  7 )  )  ; 

vise  ratio  =  flipud (parameter ( 

r 

,  8 )  )  ; 

gamma 

=  flipud (parameter ( 

r 

,10) ) ; 

s 

=  flipud (parameter ( 

r 

,11) ) ; 

mu  turb 

mu 

rho  1  = 

=  mu  lam.*visc  ratio; 
=  mu  lam  +  mu  turb; 

(  rho ( 1 : end-1 , 1 : end-1 )  + 

rho (2 : end, 1 : end-1 ) 

rho  r  = 

(  rho ( 1 : end-1 , 2 : end)  + 

rho (2 : end, 2 : end) 

)  /  2 

rho  t  = 

(  rho ( 1 : end-1 , 1 : end-1 ) 

+ 

rho ( 1 : end-1 , 2 : 

end) 

rho  b  = 

(  rho (2 : end, 1 : end-1 )  + 

rho (2 : end, 2 : end) 

)  /  2 

u _ 1  =  ( 

u ( 1 : end-1 , 1 : end-1 )  +  u(2 

end,l:end-l)  ) 

/2; 

u  r  =  ( 

u ( 1 : end-1 , 2 : end)  +  u(2 

end, 2 : end)  ) /2 ; 

u  t  =  ( 

u ( 1 : end-1 , 1 : end-1 )  +  u(l 

end-1, 2:end)  ) 

/2; 

u  b  =  ( 

u (2 : end, 1 : end-1 )  +  u(2 

end, 2 : end)  ) /2  ; 

u  c  =  ( 

ul  +  ur  +  ut  +  ub) 

/  4 ; 

v _ 1  =  ( 

v ( 1 : end-1 , 1 : end-1 )  +  v(2 

end, 1 : end-1 )  ) 

/  2  ; 

v  r  =  ( 

v ( 1 : end-1 , 2 : end)  +  v(2 

end, 2 : end)  ) /2 ; 

v  t  =  ( 

v ( 1 : end-1 , 1 : end-1 )  +  v(l 

end-1, 2:end)  ) 

/  2 ; 

v_b  =  ( 

v (2 : end, 1 : end-1 )  +  v(2 

end, 2 : end)  ) /2  ; 

v  c  =  ( 

vl+vr  +  vt  +  vb 

/  4 ; 
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T  1  =  (  T ( 1 : end-1 , 1 : end-1 )  +  T (2 : end, 1 : end-1 )  )/2; 

T  r  =  (  T ( 1 : end-1 , 2 : end)  +  T (2 : end, 2 : end)  )/2; 

T  t  =  (  T ( 1 : end-1 , 1 : end-1 )  +  T ( 1 : end-1 , 2 : end)  )/2; 

T  b  =  (  T (2 : end, 1 : end-1 )  +  T (2 : end, 2 : end)  )/2; 

T_c  =  (  T_1  +  T_r  +  T_t  +  T  Jo  ) / 4 ; 

M  1  =  (  M ( 1 : end-1 , 1 : end-1 )  +  M (2 : end, 1 : end-1 )  )/2; 

M  r  =  (  M ( 1 : end-1 , 2 : end)  +  M (2 : end, 2 : end)  )/2; 

M  t  =  (  M ( 1 : end-1 , 1 : end-1 )  +  M ( 1 : end-1 , 2 : end)  )/2; 

M  b  =  (  M (2 : end, 1 : end-1 )  +  M (2 : end, 2 : end)  )/2; 

M_c  =  (  M_1  +  M_r  +  M_t  +  M  Jo  ) / 4 ; 

gamma  1  =  (  gamma ( 1 : end-1 , 1 : end-1 )  +  gamma (2 : end, 1 : end-1 )  )/2; 

gamma  r  =  (  gamma ( 1 : end-1 , 2 : end)  +  gamma (2 : end, 2 : end)  )/2; 
gamma  t  =  (  gamma ( 1 : end-1 , 1 : end-1 )  +  gamma (1  rend- 1, 2  rend)  )/2; 
gamma  b  =  (  gamma (2 r end, 1 r end-1 )  +  gamma (2 r end, 2 r end)  )/2; 
gamma  c  =  (  gamma_l  +  gamma  r  +  gamma  t  +  gamma  Jo  ) / 4 ; 

si  =  (  s ( 1 r end-1 , 1 r end-1 )  +  s (2 r end, 1 r end-1 )  )/2; 

s_r  =  (  s ( 1 r end-1 , 2 r end)  +  s (2 r end, 2 r end)  )/ 2; 
s_t  =  (  s ( 1 r end-1 , 1 r end-1 )  +  s ( 1 r end-1 , 2 r end)  )/2; 
sjo  =  (  s (2 r end,  1 r end-1 )  +  s (2 r end, 2 r end)  )/2; 

mu  1  =  (  mu ( 1 r end-1 , 1 r end-1 )  +  mu (2 r end, 1 r end-1 )  )/2; 

mu  r  =  (  mu ( 1 r end-1 , 2 r end)  +  mu (2 r end, 2 r end)  )/2; 
mu  t  =  (  mu ( 1 r end-1 , 1 r end-1 )  +  mu ( 1 r end-1 , 2 r end)  )/2; 
mu  b  =  (  mu (2 r end, 1 r end-1 )  +  mu (2 r end, 2 r end)  )/2; 

dx  =  (  x  grid ( 1 r end-1 , 2 r end)  -  x  grid ( 1 r end-1 , 1 r end-1 )  ); 

y_t  =  (  y_gr id ( 1 r end-1 , 1 r end-1 )  +  y_grid ( 1 r end-1 , 2 r end)  )/2; 
y_b  =  (  y_grid (2 r end, 1 r end-1 )  +  y_grid (2 r end, 2 r end)  )/2; 
y_c  =  ( y_t  +  y_b ) / 2 ; 
dy  =  y_t  -  y_b; 
circ_c  =  2*pi*y__c; 


Total  Entropy  Generation 


S_l. 

S_r 

S_t_1 

S_b_l 

S_t_r 

S_b_r 

S 


zeros (n_y  nodes-l,n  x  nodes-1); 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes-l , n_x_nodes-l )  ; 
zeros (n_y_nodes-l , n_x_nodes-l )  ; 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes , n_x_nodes ) ; 


Area_l 

Area^r 

Area_t 
Area_b 
mdot_l 
mdot^r 
mdot_t 
mdot  b 


pi* (  y_grid ( 1 r end-1 , 1 r end-1 ). A2  - 
y_grid (2 r end, 1 r end-1 ) . A2  ); 
pi* (  y_grid ( 1 r end-1 , 2 r end) . A2  - 
y_grid (2 r end, 2 r end) . A2  ); 
2*pi*y_t . *dx; 

2*pi*y_b. *dx; 
rho_l . *u_l . *Area_l ; 
rho_r . *u_r . *Area_r ; 
rho  t.*v  t.*Area  t; 
rho  b.*v  b.*Area  b; 
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S_c  =  (  s_r.*mdot_r  -  s_l . *mdot_l  )  +  (  s_t.*mdot_t  - 
s_b.*mdot_b  ); 

S_l(:,2:end)  =  (  S  c (:,  1 : end-1 )  +  S_c(:,2:end)  )/2; 

S_1 ( : , 1 )  =  S_c ( : , 1 )  -  (  S_1 ( : , 2 )  -  S_c ( : , 1 )  ); 

S_r (:, 1 : end-1 )  =  S_1 ( : , 2 : end) ; 

S_r(r,end)  =  S_c(end)  +  (  S_c ( r , end)  -  S_l(:,end)  ); 

S_t_l (2 : end,  : )  =  (  S_1 (1 : end-1 ,: )  +  S_1 (2  rend,  :)  )/2; 

S_t_l ( 1 , : )  =  S_1 ( 1 , : )  +  (  S_1 ( 1 , : )  -  S_t_l(2,r)  ); 

S_b_l (1 : end-1 ,: )  =  S_t_l (2 : end, : ) ; 

S_b_l (end, :)  =  S_1 (end, : )  -  (  S_b_l (end-1 ,: )  -  S_1 (end, : )  ); 

S_t_r (2 : end, : )  =  (  S_r (1 : end-1 ,: )  +  S_r(2:end, :)  )/2; 

S_t_r ( 1 , r )  =  S_r ( 1 , : )  +  (  S_r(l,r)  -  S_t_r(2,r)  ); 

S_b_r (end, end)  =  S_r(end,end)  -  (  S_b_r (end-1 , end)  - 
S_r (end, end)  ) ; 

S (1  rend- 1, 1  rend- 1)  =  S_t_l; 

S (end, 1 r end-1 )  =  S_b_l (end, r ) ; 

S ( 1 r end-1 , end)  =  S_t_r ( : , end) ; 

S (end, end)  =  S_b_r (end, end) ; 

S  ( : ,  1 )  =  S  ( : ,  2 )  -  (  S  ( r ,  3)  -  S  ( r ,  2 )  )  ; 

S(:,end)  =  S(r, end-1)  +  (  S(r, end-1)  -  S(r, end-2)  ); 


parameter (:,:, 12 )  =  flipud(S); 


Entropy  Gen.  Due  to  Friction 


tau_x_t 

tau_x_b 

tau_y_l 

tau_y_r 


zeros (n_y  nodes-l,n  x  nodes-1); 
zeros (n_y_nodes-l , n_x_nodes-l ) ; 
zeros (n_y_nodes-l , n_x_nodes-l ) ; 
zeros (n_y  nodes-1, n  x  nodes-1); 


S^f^! 

S_f_r 

S_f_t_! 

S_f_b_l 

S_f_t_r 

S_f_b_r 

S  f 


zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes-l , n_x_nodes-l ) ; 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes-l , n_x_nodes-l )  ; 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes , n_x_nodes )  ; 


Second  Order  Accurate 


tau_x_t ( 1 , r )  =  mu_t ( 1 , r ) . * (  u_t ( 1 , : )  -  u_b ( 1 , r )  ) . /dy ( 1 , r ) ; 

tau_x_t (2  r  end,  r )  =  mu_t (2  r  end,  : )  . * (  u_t ( 1 r  end-1 ,  : )  - 

u_b (2  rend,  : )  )  . / (  dy (1 r end-1,  : )  +  dy (2  rend,  : )  ) ; 
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tau_x_b ( 1 : end- 1 ,  : )  =  mu_b ( 1 : end- 1 , : ) .  *  (  u_t ( 1 : end- 1 ,  : )  - 
u_b(2:end, :)  )./(  dy ( 1 : end-1 ,: )  +  dy(2:end, :)  ); 

tau_x_b (end, : )  =  mu_b (end, : ) . * (  u_t (end, : )  -  u_b (end, : ) 

) . /dy (end, : ) ; 

tau_y_l  (bl)  =  mu_l  (  :  ,  1 )  .  *  (  v_r  (  :  ,  1 )  -  v_l  (  :  ,  1 )  )  .  /dx  (  :  ,  1 )  ; 
tau_y_l ( : , 2 : end)  =  mu_l ( : , 2 : end) . * (  v_r ( : , 2 : end)  - 

v_l (:, 1 : end-1 )  )./ (dx (:, 1 : end-1 )  +  dx(:,2:end)  ); 

tau_y_r ( : , 1 : end-1 )  =  mu_r ( : , 1 : end-1 ) . * (  v_r ( : , 2 : end)  - 

v_l (:, 1 : end-1 )  )./(  dx (:, 1 : end-1 )  +  dx(:,2:end)  ); 

tau_y_r ( : , end)  =  mu_r ( : , end) . * (  v_r ( : , end)  -  v_l ( : , end) 

) . /dx ( : , end) ; 


S_f_c  =  circ_c.*(  l./T_c.*(  dx . *tau_x_t . * (u_t  -  u_c)  +  ... 

dx. *tau  x  b.*(u  c  -  u  b)  +  dy.*tau  y_r.*(v  r  -  v_c)  + 

dy .  *tau_y_l . * (v_c  -  v_l)  )  ); 

S_f_l ( : , 2 : end)  =  (  S_f_c (:, 1 : end-1 )  +  S_f_c ( : , 2 : end)  )/2; 
S_f_l ( : , 1 )  =  S_f_c ( : , 1 )  -  (  S_f_l ( : , 2 )  -  S_f_c ( : , 1 )  ) ; 

S_f_r ( : , 1 : end-1 )  =  S_f_l ( : , 2 : end) ; 

S_f_r(:,end)  =  S_f_c(end)  +  (  S_f_c(:,end)  -  S_f_l ( : , end)  ); 

S_f_t_l (2 : end, : )  =  (  S_f_l (1 : end-1, : )  +  S_f_l (2 : end, : )  )/2; 
S_f_t_l ( 1 , : )  =  S_f_l ( 1 , : )  +  (  S_f_l ( 1 , : )  -  S_f_t_l(2,:)  ); 

S_f_b_l ( 1 : end-1 ,: )  =  S_f_t_l(2:end,:); 

S_f_b_l (end, :)  =  S_f_l (end, :)  -  (  S_f_b_l (end-1 , :)  - 
S_f_l (end, : )  ); 

S_f_t_r (2 : end, : )  =  (  S_f_r (1 : end-1 ,: )  +  S_f_r (2 : end, : )  )/2; 

S_f_t_r ( 1 , : )  =  S_f_r ( 1 , : )  +  (  S_f_r(l,:)  -  S_f_t_r(2,:)  ); 

S_f_b_r (end, end)  =  S_f_r (end, end)  -  (  S_f_b_r (end-1 , end)  - 

S_f_r (end, end)  ); 

S_f (1 : end-1 , 1 : end-1 )  =  S_f  t  1; 

S_f (end, 1 : end-1 )  =  S_f_b_l(end,:); 

S_f ( 1 : end-1 , end)  =  S_f_t_r ( : , end) ; 

S_f(end,end)  =  S_f_b_r (end, end) ; 

S_f ( : , 1 )  =  S_f(:,2)  -  (  S_f(:,3)  -  S_f(:,2)  ); 

S_f(:,end)  =  S_f ( : , end-1)  +  (  S _ f ( : , end-1)  -  S_f ( : , end-2)  ); 

parameter (:,:, 13)  =  flipud(Sf); 
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Entropy  Gen.  Due  to  conv.  HT 


conv  HT_1  =  zeros (n  y_nodes-l,n  x  nodes-1); 
conv  HT  r  =  zeros (n  y_nodes-l,n  x  nodes-1); 
conv  HT  t  =  zeros (n  y_nodes-l,n  x  nodes-1); 
conv  HT  b  =  zeros (n  y_nodes-l,n  x  nodes-1); 


S_conv_HT  1 
S_conv_HT_r 
S_conv_HT  t_l 
S_conv_HT  b_l 
S_conv_HT_t_r 
S_conv_HT  b  r 
S  conv  HT 


zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes-l , n_x_nodes-l )  ; 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes-l , n_x_nodes-l )  ; 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes , n_x_nodes )  ; 


Pr_lam  =  0.72; 
Pr  turb  =  0.90; 


k  lam  =  (  gamma. *P.*mu  lam  )./(  (gamma  -  1 ) . *rho . *T*Pr  lam  ); 
k_turb  =  (  gamma. *P.*mu  turb  ) ./(  (gamma  -  1) . *rho . *T*Pr_turb  ); 
k  =  k  lam  +  k  turb; 


k  1  =  (  k ( 1 : end-1 , 1 : end-1 )  +  k (2 : end, 1 : end-1 )  )/2; 
k  r  =  (  k ( 1 : end-1 , 2 : end)  +  k (2 : end, 2 : end)  )/2; 
k  t  =  (  k (1 :end-l, 1 :end-l)  +  k ( 1 : end-1 , 2 : end)  )/2; 
k  b  =  (  k (2 : end, 1 : end-1 )  +  k (2 : end, 2 : end)  )/2; 
conv  HT_r ( : , 1 : end-1 )  =  k_r ( : , 1 : end-1 )  . * (  T  r ( :  ,  2 : end)  - 

T_1 (:, 1 : end-1)  )./(  dx (:, 1 : end-1 )  +  dx(:,2:end) 
conv  HT_r ( : , end)  =  k  r ( : , end) . * (  T  r ( : , end)  -  T  1 ( : , end-2 ) 

dx ( : ,  end-2)  +  dx ( : , end-1)  +  dx ( : , end)  ) 


)  ; 

)  •/  ( 


conv_HT_l ( : , 1 )  =  k_l(:,l).*(  T_r ( : , 1 )  -  T_1 ( : , 1 )  )./(  dx ( : , 1 )  ); 

conv  HT_1 ( : , 2 : end)  =  k_l ( : , 2 : end) . * (  T  r ( : , 2 : end)  -  T_1 ( : , 1 : end- 

1)  )./(  dx (:, 1 : end-1 )  +  dx(:,2:end)  ); 


conv_HT_t (1, : )  =  k_t(l,:).*(  T_t ( 1 , : )  -  T_l(l,:)  )/(  dy(l,:)  ); 

conv_HT_t (2 : end, : )  =  k_t (2 : end, : ) . * (  T_t ( 1 : end-1 ,: )  - 

T_b(2:end, :)  )./(  dy ( 1 : end-1 ,: )  +  dy(2:end, :)  ); 

conv  HT  b (1 :end-l,  : )  =  k  b (1 :end-l,  : )  . * (  T  t (1 :end-l,  : )  - 

T_b(2:end, :)  )./(  dy ( 1 : end-1 ,: )  +  dy(2:end, :)  ); 

conv_HT_b (end,  : )  =  k_b (end,  : )  . * (  T_t (end,  : )  -  T_b (end,  : )  )  . / ( 

dy (end, : )  ) ; 

S_conv_HT_c  =  circ_c.*(  l./T_c.*(  dy.*conv_HT  r  -  dy.*conv  HT_1 
+  dx.*conv  HT  t  -  dx.*conv  HT  b  )  ); 


S_conv_HT  l(:,2:end)  =  (  S_conv  HT_c (:, 1 : end-1 )  + 

S_conv_HT_c ( : , 2 : end)  ) /2  ; 

S_conv_HT  1 ( :  ,  1 )  =  S_conv  HT  c(:,l)  -  (  S  conv  HT_1(:,2)  - 

S_conv_HT_c ( : , 1 )  )  ; 

S  conv  HT  r(:,l: end-1)  =  S  conv  HT  l(:,2:end); 
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S_conv_HT  r ( :  ,  end)  =  S_conv  HT_c(end)  +  (  S_conv_HT  c ( : , end)  - 

S_conv_HT_l ( : , end)  )  ; 

S  conv  HT  t_l(2:end, :)  =  (  S  conv  HT_1 (1 : end-1 ,: )  + 

S_conv_HT_l (2 : end, : )  )/2; 

S_conv  HT  t_l(l,:)  =  S_conv  HT_1(1,:)  +  (  S_conv_HT  1(1,:)  - 

S_conv_HT_t_l (2, : )  ) ; 

S_conv_HT  b_l ( 1 : end-1 ,: )  =  S  conv  HT  t_l (2 : end, : ) ; 

S_conv  HT  b_l  (end,  :  )  =  S__conv_HT  1  (end,  :  )  -  ( 

S_conv  HT  b  1 (end-1, : )  -  S  conv  HT  1 (end, : )  ) ; 

S_conv_HT  t  r(2:end, :)  =  (  S_conv  HT_r (1 : end-1 ,: )  + 

S_conv_HT_r (2 : end,  : )  )/2; 

S_conv  HT  t_r(l, :)  =  S_conv  HT_r(l,:)  +  (  S_conv_HT  r(l,:)  - 

S_conv_HT_t_r (2, : )  ) ; 

S_conv_HT  b_r(end,end)  =  S_conv  HT_r (end, end)  -  ( 

S_conv_HT  b  r (end-1,  end)  -  S_conv^HT  r(end,end)  ); 

S_conv  HT ( 1 : end-1 , 1 : end-1 )  =  S_conv  HT  t_l; 

S_conv  HT (end,  1 : end-1 )  =  S_conv_HT  b  1 (end,  : ) ; 

S_conv  HT ( 1 : end-1 , end)  =  S_conv  HT  t_r(:,end); 

S_conv  HT (end, end)  =  S_conv  HT  b  r (end, end) ; 

S_conv_HT ( : , 1 )  =  S_conv  HT(:,2)  -  (  S  conv  HT(:,3)  - 
S_conv_HT ( : , 2 )  ); 

S_conv_HT ( : , end)  =  S_conv_HT (:, end-1 )  +  (  S_conv_HT ( : , end-1 )  - 

S_conv_HT ( : , end-2 )  ) ; 

parameter (:,:, 14 )  =  flipud (S_conv_HT) ; 

% - Entropy  Gen.  Due  to  non-ideal  HT - 


S_nideal  HT  =  zeros (n_y  nodes, n  x  nodes); 
parameter (:,:, 15)  =  flipud (S_nideal_HT) ; 


if  b  ==  19  ||  b  ==  20 
S_nideal_HT_l  = 
S^nideal  HT  r  = 
S^nideal^HT  t_l  = 
S^nideal  HT  b_l  = 
S^nideal  HT  t  r  = 
S^nideal  HT  b_r  = 
S  nideal  HT  = 


zeros (n_y_nodes-l , n_x_nodes-l ) ; 
zeros (n_y  nodes-l,n  x  nodes-1); 
zeros (n_y_nodes-l , n_x_nodes-l ) ; 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y  nodes-1, n  x  nodes-1); 
zeros (n_y_nodes-l , n_x_nodes-l )  ; 
zeros (n_y_nodes , n_x_nodes ) ; 


T  t_c  =  T  c . * (  1+ (gamma  c  -  1)./2.*M  c .  A2  ); 


V  =  pi*r_iA2* (l/2*L_b) ; 

dV  =  pi* (  y_t . A 2  -  y_b.A2  ).*dx; 


S  nideal  HT  c 


(l/2*Qdot) * (dV./V) .* (l./T  c 


l./T  t  c  ) ; 
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S  nideal_HT_l ( : , 2 : end)  =  (  S  nideal_HT_c (:, 1 : end-1 )  + 

S_nideal_HT_c ( :  ,  2 : end)  )/2; 
S_nideal_HT_l ( : , 1 )  =  S_nideal_HT_c ( :  ,  1 )  -  ( 

S  nideal_HT  1(:,2)  -  S_nideal  HT_c(:,l)  ); 

S  nideal  HT  r(:,l: end-1)  =  S  nideal_HT  l(:,2:end); 

S  nideal_HT_r ( : , end)  =  S  nideal_HT_c ( : , end)  +  ( 

S_nideal_HT_c ( : , end)  -  S_nideal_HT_l ( : , end)  ) ; 

S  nideal  HT  t_l(2:end, :)  =  (  S_nideal  HT_1 (1 : end-1 ,: )  + 

S_nideal_HT_l (2 : end,  : )  )/2; 

S  nideal_HT  t_l(l,:)  =  S  nideal_HT_l ( 1 , : )  +  ( 

S  nideal_HT  1(1,:)  -  S_nideal  HT  t  1(2,:)  ); 

S  nideal  HT  b_l (1 :end-l, : )  =  S_nideal  HT  t  1(2: end, :); 

S  nideal_HT  b_l (end, :)  =  S_nideal  HT_1 (end, :)  -  ( 
S_nideal  HT  b  1 (end-1,:)  -  S_nideal  HT_1 (end, :)  ); 

S  nideal  HT  t  r(2:end, :)  =  (  S_nideal  HT_r (1 : end-1 ,: )  + 

S_nideal_HT_r (2 : end,  : )  )/2; 

S  nideal^HT  t_r(l,:)  =  S  nideal_HT_r ( 1 , : )  +  ( 

S  nideal_HT  r(l,:)  -  S_nideal  HT  t  r(2,:)  ); 

S  nideal  HT  b  r (end, end)  =  S  nideal_HT  r (end, end)  -  ( 
S_nideal  HT^t  r (end, end)  -  S_nideal  HT_r (end, end)  ); 

S  nideal_HT ( 1 : end-1 , 1 : end-1 )  =  S_nideal_HT  t_l; 

S  nideal  HT (end,  1 : end-1 )  =  S  nideal_HT  b_l (end,  :); 

S  nideal^HT ( 1 : end-1 , end)  =  S  nideal_HT  t_r ( : , end) ; 

S  nideal  HT (end, end)  =  S  nideal_HT  b_r (end, end) ; 

parameter (:,:, 15)  =  flipud (S_nideal_HT) ; 

end 


% - Entropy  Gen.  Due  to  ideal  HT - 

S_ideal_HT  =  zeros (n_y  nodes, n  x  nodes); 
parameter (:,:, 1 6 )  =  flipud (S_ideal_HT) ; 
if  b  ==  19  ||  b  ==  20 

S  ideal_HT  1  =  zeros (n  y_nodes-l,n  x  nodes-1); 

S  ideal_HT  r  =  zeros (n  y_nodes-l , n_x  nodes-1); 

S  ideal_HT  t_l  =  zeros (n  y__nodes-l,n  x  nodes-1); 

S  ideal_HT  b_l  =  zeros (n  y  nodes-1, n_x  nodes-1); 

S  ideal_HT  t_r  =  zeros (n  y_nodes-l , n_x  nodes-1); 

S  ideal_HT  b  r  =  zeros (n  y_nodes-l,n  x  nodes-1); 

S  ideal_HT  =  zeros (n  y_nodes,n  x  nodes); 

S_ideal_HT_c  =  (l/2*Qdot) * (dV. /V) . /T_t_c; 

S  ideal_HT  l(:,2:end)  =  (S_ideal  HT_c (:, 1 : end-1 )  + 

S_ideal_HT_c ( : ,  2 : end)  )/2; 
S_ideal_HT_l ( : , 1 )  =  S_ideal_HT_c ( : , 1 )  -  ( 

S  ideal  HT  1(:,2)  -  S  ideal  HT  c(:,l)  ); 
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S  ideal_HT  r(:,l: end-1)  =  S_ideal  HT_1  (  :  ,  2  :  end)  ; 

S  ideal_HT  r ( : , end)  =  S  ideal_HT  c ( : , end)  +  ( 

S_ideal_HT_c ( : , end)  -  S_ideal_HT_l ( : , end)  ) ; 

S  ideal_HT  t_l(2:end, :)  =  (  S_ideal_HT  l(l:end-l,:)  + 

S_ideal_HT_l (2 : end,  : )  )/2; 

S_ideal_HT_t_l (1, : )  =  S_ideal_HT_l (1, : )  +  ( 

S_ideal_HT_l (1, : )  -  S_ideal_HT_t_l (2 , : )  ); 

S  ideal_HT  b__l  (1 :  end-1 ,: )  =  S_ideal_HT  t_l  (2  :  end,  : )  ; 

S  ideal_HT  b__l  (end,  :)  =  S_ideal_HT  1  (end,  :)  -  ( 

S_ideal_HT_b_l (end-1, : )  -  S_ideal_HT_l (end, : )  ) ; 

S  ideal_HT  t_r(2:end, :)  =  (  S_ideal_HT  r(l:end-l, :)  + 

S_ideal_HT_r (2 : end, : )  )/2; 

S_ideal_HT_t_r ( 1 , : )  =  S_ideal_HT_r ( 1 , : )  +  ( 

S_ideal_HT_r ( 1 , : )  -  S_ideal_HT_t_r (2 , : )  ); 

S  ideal_HT  b  r (end, end)  =  S_ideal  HT^r (end, end)  -  ( 
S_ideal_HT_t_r (end, end)  -  S_ideal_HT_r (end, end)  ); 

S  ideal_HT ( 1 : end-1 , 1 : end-1 )  =  S_ideal  HT__t  1; 

S  ideal_HT (end, 1 : end-1 )  =  S  ideal  HT  b  1 (end, : ) ; 

S  ideal_HT ( 1 : end-1 , end)  =  S  ideal  HT  t  r ( : , end) ; 
S_ideal_HT (end, end)  =  S_ideal_HT_b_r (end, end) ; 

parameter (:,:, 1 6 )  =  flipud (S_ideal_HT) ; 

end 


% - Entropy  Gen.  Due  to  HT - 

S_HT  =  S_conv_HT  +  S_nideal_HT  +  S_ideal_HT; 
parameter (:,:, 17 )  =  flipud (S_HT) ; 

% - Entropy  Gen.  Due  to  Shocks - 

S_shock  =  S  -  (S_f  +  S_HT); 
parameter (:,:, 1 8 )  =  f lipud ( S_shock) ; 
for  p  =  1 : nn_parameters 

dlmwrite (write_label , parameter ( : , : , p) , ' -append ' , 

'delimiter',  1  ;,  'precision',  '%.15f') 

end 

end 


end 


o\°  o\°  o\° 
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% - Name  File - 

for  p  =  1 : nn_parameters 
if  p  ==  1 

dlmwrite ( ' timhis  solution  name . nam ' , parameter^name (p, : ) , ' ' ) 

else 

dlmwrite (' timhis  solution  name . nam pa ramete rename (p, - 
append',  'delimiter',  '') 

end 

end 

fprintf (' timhis  solution  name. nam  File  Created  \n') 
toe 

fprintf (' Creation  of  timhis  X  solution. f  Files  Completed  \n  \n') 

end 


timhis  animation  File 


if  Create  Contour  Animation  ==  1 
Contour  =  6; 
for  t  =  l:n  timhis 

read_label  =  [ ' timhis_ ' , int2str ( t ) , ' _solution . f ' ] ; 
write  label  =  'timhis  animation . f ' ; 

fprintf (' Start  Creation  of  timhis  solution  animation:  %s 
\n ' , read_label ) 

timhis  solution  =  dlmread(read  label); 
if  t  ==  1 

dlmwrite (write_label , n_timhis ,  '  ') 

n  t  x  nodes  =  n  x  nodes*n  blocks  -  n  blocks  +  1; 
for  j  =  l:n  timhis 

dlmwrite (write_label ,  [n  t  x  nodes  n  y  nodes],  '-append', 
' delimiter ' ,  '  ' ) 

dlmwrite (write_label , 1 , '-append1,  'delimiter',  '  ') 

end 

end 

timhis  animate  =  zeros (n  y_nodes,n  blocks*n  x  nodes- 

n  x  nodes+1 ) ; 

for  b  =  l:n  blocks 
if  b  ==  1 

first_column  =  1; 

last_column  =  first  column  +  (n  x  nodes-1); 
first_row  =  2*n  blocks  +  2  +  (Contour-1 ) *n_y  nodes; 
last_row  =  first_row  +  n  y_nodes  -  1; 

else 

first_column  =  last_column; 

last_column  =  first  column  +  (n  x  nodes-1); 
first_row  =  last  row  +  (nn_parameters-l) *n_y  nodes  +  1; 
last_row  =  first_row  +  n  y_nodes  -  1; 

end 

timhis  animate  (:,  first_column  :  last__column)  = 

timhis  solution ( first  row: last  row, : ) ; 


end 
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end 


dlmwrite (write  label, timhis  animate, ' -append ' , 
',  'precision',  '%.8e') 


' delimiter ' , 


% - Grid  Construction  Per  Block - 

X  block  =  zeros (n_y  nodes, n  x  nodes, n  blocks); 

Y  block  =  zeros (n_y  nodes, n  x  nodes, n  blocks); 

for  b  =  l:n  blocks 

fprintf (' Begin  Grid  Construction  for  Block:  %4.0f  \n',b) 
if  b  ==  1 

first  column  =  1; 

last_column  =  first_column  +  (n  x  nodes-1); 
first  row  x  =  n  blocks  +  2; 

last_row_x  =  first  row  x  +  (n_y_nodes-l ) ; 

first  row_y  =  last_row  x  +  1; 

last_row_y  =  first  row_y  +  (n_y_nodes-l ) ; 

else 

first^column  =  last  column; 

last_column  =  first_column  +  (n  x  nodes-1); 

first  row  x  =  last_row_y  +  1; 

last^row^x  =  first  row_x  +  (n_y_nodes-l ) ; 

first  row_y  =  last_row  x  +  1; 

last_row_y  =  first  row_y  +  (n_y  nodes-1); 

end 

X ( : , f irst^column : last_column)  = 

0 . 25*grid (first  row  x:last_row_x, :) ; 

Y  (:,  first  column  :  last__column)  = 

0 . 25*grid (first  row^yrlast  row  y, : ) ; 

end 

XY  =  [  X ; Y ] ; 

dlmwrite (' grid  animate . grd ',  n  timhis,  '  ') 

for  t  =  l:n  timhis 

dlmwrite (' grid  animate . grd ', [n  t  x  nodes,  n  y  nodes ],' -append  , 

' delimiter ' ,  '  ' ) 

end 

for  t  =  l:n  timhis 

dlmwrite ( 'grid  animate . grd ' , XY,  ' -append ' ,  ' delimiter '  ,  ' 

',  'precision',  ' % . 8  f ' ) 

end 

fprintf (' Animate  Grid  File  Complete  \n') 

dlmwrite ( ' timhis  animate  name . nam 1 , parameter^name (Contour, :),'') 

fprintf (' Animate  Name  File  Complete  \n') 

toe 

fprintf (' Creation  of  timhis  animation. f  File  Completed  \n  \n') 


end 


o\°  o\°  o\° 
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Center-Line  Plot 


if  Create  CL  Plot  ==  1 

plot_parameter  =  6; 

plot_timhis_solution  =  100; 

fprintf (' Start  Creation  of  Center-Line  Plot  \n') 
timhis  solution  = 

dlmread ( [ ' timhis  ' , int2str (plot  timhis^solution) , '  solution . f ']) ; 
% - Solution - 


for  b  =  l:n  blocks 

fprintf (' Begin  Reconstruction  for  Block:  %4.0f  \n',b) 
if  b  ==  1 

first  column  =  1; 

last_column  =  first_column  +  (n  x  nodes-1); 

first  row  =  2*n  blocks  +  2  +  (plot_parameter-l ) *n_y  nodes; 

last^row  =  first  row  +  n  y  nodes  -  1; 

else 


first  column  =  last  column; 

last_column  =  first_column  +  (n  x  nodes-1); 

first  row  =  last_row  +  (nn_parameters-l) *n  y_nodes  +  1; 

last_row  =  first  row  +  n  y  nodes  -  1; 

end 


end 


CL_plot ( 1 , f irst_column : last  column)  = 

timhis  solution ( first  row. 


Grid 


CL  X  =  zeros (l,n  blocks*n  x  nodes-n  x  nodes+1); 


for  b  =  l:n  blocks 
fprintf ( ' Begin 
if  b  ==  1 

first  column 
last_column 
first  row  x 
last_row_x 
first  row_y 
last_row_y 

else 

first  column 
last_column 
first  row  x 
last_row  x 
first_row_y 
last_row_y 

end 


Construction  for 

:  1; 

:  first_column  +  (n 
n  blocks  +  2; 
first  row  x  +  (n_y 
last  row  x  +  1; 
first  row  y  +  (n_y 

:  last  column; 

:  first_column  +  (n 
last  row_y  +  1; 
first  row  x  +  (n_y 
last  row  x  +  1; 
first  row  y  +  (n_y 


Block:  %4 . Of  \n',b) 

x_nodes-l ) ; 
nodes-1 ) ; 
nodes-1 ) ; 

x_nodes-l ) ; 
nodes-1 ) ; 
nodes-1 ) ; 


Grid 


CL  X ( 1 , f irst_column : last_column)  =  0 . 25*grid (first  row  x. 


End 


o\°  o\°  o\° 
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Plot 


font_size  =  20; 

Title  =  ['Center  Line  , parameter  name (plot_parameter , :), 'vs. 

Axial  Location']; 

Y  label  =  ['Center  Line  ', parameter  name (plot_parameter,  : ) ]  ; 
figure  ( 1 ) 

plot (CL  X, CL_plot, ' LineWidth ' , 2 ) 

title (Title, ' FontName ',' Times  New  Roman ',' Fontsize ', font  size) 
xlabel (' Axial  Location  [m] ',' FontName ',' Times  New 
Roman ' , ' Fontsize ' , font  size) 

ylabel (Y  label ,' FontName ',' Times  New  Roman ',' Fontsize ', font  size) 

set (gca, ' FontName ',' Times  New  Roman ',' Fontsize ', font  size) 

toe 

fprintf (' Creation  of  Center-Line  Plot  Completed  \n  \n') 

end 


Line  Graph  Movie/Entropy  Results 


if  Create  LG_movie  ==  1 

fprintf (' Start  Creation  of  Line  Graph  Movie/Entropy  Results  \n') 
% - Asemble  Grid  into  Matrix - 

X  =  zeros (n  y  nodes,  n  blocks*n  x  nodes-n  x  nodes+1); 

Y  =  zeros (n_y  nodes, n  blocks*n  x  nodes-n_x  nodes+1); 

for  b  =  l:n  blocks 

fprintf (' Begin  Grid  Construction  for  Block:  %4.0f  \n',b) 
if  b  ==  1 

first_column  =  1; 

last_column  =  first  column  +  (n  x  nodes-1); 
first_row_x  =  n  blocks  +  2; 

last  row  x  =  first  row  x  +  (n  y  nodes-1); 

first_row_y  =  last_row  x  +  1; 

last  row_y  =  first  row_y  +  (n  y  nodes-1); 

else 

first_column  =  last  column; 

last_column  =  first  column  +  (n  x  nodes-1); 

first_row_x  =  last_row_y  +  1; 

last_row  x  =  first  row  x  +  (n  y  nodes-1); 

first  row_y  =  last^row  x  +  1; 

last_row__y  =  first  row_y  +  (n_y_nodes-l )  ; 

end 

X ( : , f irst_column : last_column)  = 

0 . 25*flipud (grid (first  row  x:last  row  x,:)); 
Y ( : , f irst_column : last_column)  = 

0 . 25*flipud (grid ( f irst_row_y : last_row_y, : ) ) ; 


end 
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Area  =  pi* (  Y ( 1 : end-1 A2  -  Y (2 : end, : ) . A2  ); 
% -  Assemble  Data  into  Matricies 


timhis_parameter  =  zeros (n_y  nodes, n  blocks*n  x  nodes- 

n  x  nodes+1 , nn_parameters ) ; 
thrust^ss  =  zeros ( 1 , n_timhis ) ; 
thrust  t  =  zeros (l,n  timhis); 
time  =  zeros (l,n  timhis); 


for  t  =  l:99:n  timhis 
timhis_solution  = 

dlmread ( [ ' timhis  ' , int2str (t) , '  solution . f ' ] ) ; 

if  t  ==  1 


first  column  =  zeros (l,n  blocks); 
last  column  =  zeros (l,n  blocks); 
first  row  =  zeros (nn_parameters , n  blocks); 
last  row  =  zeros (nn_parameters , n  blocks); 
for  b  =  l:n  blocks 
if  b  ==  1 


first  column (b)  =  1; 

last_column (b)  =  first  column (b)  +  n  x  nodes  -  1; 

else 


first  column (b)  =  last_column (b-1 ) ; 

last_column (b)  =  first  column (b)  +  n  x  nodes  -  1; 

end 

for  p  =  1 : nn_parameters 

first  row(p,b)  =  2*n  blocks  +  2  +  (p  -  l)*n  y^nodes 
+  (b  -  1 ) *nn_parameters*n  y_nodes; 
last_row (p, b)  =  f irst^row (p, b)  +  n_y_nodes  -  1; 


end 


end 


end 

for  b  =  l:n  blocks 

for  p  =  1 : nn_parameters 

timhis_parameter ( 1 : n_y  nodes , first_column (b) : last_column (b) , p)  = 
flipud (timhis  solution ( first  row (p, b) : last  row(p,b),l:n  x  nodes)); 
end 


end 


Average  of  nodal  values 


timhis_parameter_avg  =  (timhis_parameter (1 :end-l, : , : )  + 

timhis_parameter (2 : end,  :  ,  : ) ) / 2 ; 


rho_avg 

u_avg 

P_avg 

mu_lam 

vise  ratio 

s_avg 

S_f_avg 

S_conv  HT  avg 

S_nideal_HT  avg 

S  ideal  HT  avg 


timhis  parameter  avg ( 

r 

,D  ; 

timhis  parameter  avg ( 

r 

,2)  ; 

timhis  parameter  avg ( 

r 

,4)  ; 

timhis  parameter  avg ( 

r 

,7)  ; 

timhis  parameter  avg ( 

r 

,8)  ; 

timhis  parameter  avg ( 

r 

,11) 

timhis  parameter  avg ( 

r 

,13)  ; 

timhis  parameter  avg ( 

r 

,14) 

timhis  parameter  avg ( 

r 

,15) 

timhis  parameter  avg ( 

r 

,16)  ; 

,17)  ; 


S_HT_avg 
mu  turb 
mu 


=  timhis_parameter  avg ( : , : 
=  mu  lam.*visc  ratio; 

=  mu  lam  +  mu  turb; 


S_f_plot 
S_conv_HT_plot 
S_nideal  HT_plot 
S_ideal_HT_plot 
S_HT_plot 


sum ( S_f_avg) ; 
sum ( S_conv^HT  avg); 
sum(S  nideal_HT  avg); 
sum(S  ideal_HT  avg); 
sum ( S_HT_avg) ; 


Axial  integration  of  entropy  generation 


last  x  node  =  n  blocks*n  x  nodes  -  n  blocks  +  1; 


for  i  =  2: last  x  node 
S_f_plot ( i ) 
S_conv_HT_plot (i) 

S  nideal_HT_plot ( i ) 

S_ideal_HT_plot (i) 

S_HT_plot (i) 

end 


S_f_plot (i-1 )  +  S_f_plot(i); 
S_conv_HT_plot ( i-1 )  + 
S_conv_HT_plot (i) ; 

S_nideal  HT_plot(i-l)  + 
S_nideal  HT_plot(i); 
S_ideal_HT_plot (i-1)  + 
S_ideal_HT_plot (i) ; 

S_HT_plot (i-1)  +  S_HT_plot ( i ) 


mdot  =  rho  avg.*u  avg.*Area; 

S_plot  =  sum (mdot . *s_avg) ; 

S_plot(l)  =  S_plot(2); 

S_plot  =  S_plot  -  S_plot(l); 

S_plot  =  S_plot  -  S_plot(501)  +  (S__f_plot  (501)  + 
S_HT_plot (501) )  ; 

S_shock_plot  =  S_plot  -  (S_f_plot  +  S  conv  HT_plot  + 
S_nideal  HT_plot  +  S^ideal  HT_plot) ; 

axial  mdot  =  sum(mdot); 

axial  Sdot  avg  =  sum (mdot . *s_avg) ; 

axial  Hdot  t_avg  =  sum (mdot . *h_t  avg) ; 

xaxis  =  X(end, :); 

y  axis_l  =  axial_mdot/axial  mdot(l); 
y  axis_2  =  axial_Sdot  avg; 

y  axis_3  =  axial_Hdot_t  avg/max (axial  Hdot  t_avg) ; 


movie  matrix 


[x  axis',y  axis  1']; 
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Thrust  Calculations 


if  t  ==  1 

y_l_avg_x  =  1/2* (  Y(l,l:end-1)  +  Y(l,2:end)  ); 
y_2_avg_x  =  1/2* (  Y(2,l:end-1)  +  Y (2, 2  rend)  ); 
x  1  =  X ( 1 , 1 : end-1 ) ; 
x  2  =  X ( 1 , 2 : end) ; 
y_l  =  Y  ( 1 , 1 : end-1 ) ; 
y_2  =  Y ( 1 , 2 : end) ; 

dl  =  sqrt (  (x_2  -  x_l) .  A2  +  (y_2  -  y_l) .  A2  ); 

=  atan  (  (y_2  -  y_l)  ./ (x_2  -  x_l)  ) ; 


=  1/2* (  timhis_parameter ( 1 , 1 : end-1 , 4 )  + 

timhis_parameter ( 1 , 2 : end,  4  )  ); 

=  1/2*  (  timhis_parameter ( 1 , 1 : end-1 , 2 )  + 

timhis_parameter ( 1 , 2 : end,  2 )  ); 

=  1/2*  (  timhis_parameter (2 , 1 : end-1 , 2 )  + 

timhis_parameter (2 , 2 : end,  2 )  ); 

=  l/2*(  mu ( 1 , 1 : end-1 )  +  mu(l,2:end)  ); 

=  (u  2  avg  x  -  u_l  avg  x) / (y_l  avg  x  -  y  2^avg  x) ; 

=  mu_avg_x*du_dy; 
thrust_ss (t)  =  sum ( (mdot (:, end) . *u_avg (:, end)  + 

P_avg ( : , end) . *Area ( : , end) )  - 
(mdot ( : , 1 ) . *u_avg ( : , 1 )  + 

P_avg ( : , 1 ) . *Area ( : ,  1 )  )  )  ; 

thrust_t  P  =  sum(  2*pi*P  avg  x.*dl.*y  l_avg  x . *sin (alpha)  ); 
thrust  t  tau  =  sum(  -pi*tau . * (x^2  -  x  l).*(y_l  + 

y_2 ). *cos (alpha)  ); 

thrust_t(t)  =  thrust_t_P  +  thrust_t_tau; 
time(t)  =  1.0e-5*t; 

% -  Begin  Construction  of  Movie  File  for  Techplot  360  - 

line_l  =  ' TITLE="Title" ' ; 

line_2  =  ' VARIABLES="X"  "Y" '  ; 

line_3  =  'ZONE  1=2400  F=POINT'; 

line_end  =  ['TEXT  X=10  Y=90  T="Time  = 

' , num2str ( (t) * . 00001, '%4.5f'),  ,  ZN=  ' , int2str (t) ] ; 

if  t  ==  1 

dlmwrite ( ' Line  Graph  Movie . txt ' , line_l ,  '  '  ) 
dlmwrite ( ' Line  Graph_Movie . txt ', line  2,'- 
append' , 'delimiter' , ' ' ) 
dlmwrite (' Line  Graph_Movie . txt ', line  3,'- 
append' , 'delimiter' , ' ' ) 

dlmwrite (' Line  Graph  Movie . txt ', movie  matrix, '- 

append ' ,  ' delimiter ', ' \t ' ,  'precision',  ' % . 8  f ' ) 
dlmwrite (' Line  Graph  Movie . txt 1 , line  end, ' - 
append ; , ; delimiter ' , ' ' ) 


alpha 

end 

P  avg  x 

u  l_avg  x 

u  2  avg  x 

mu  avg  x 

du_dy 

tau 


else 
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dlmwrite ( ' Line_Graph  Movie . txt ' , line_l , ' - 
append' , 'delimiter' , ' ' ) 
dlmwrite ( ' Line_Graph  Movie . txt ' , line_2 , ' - 
append' , 'delimiter' , ' ' ) 
dlmwrite ( ' Line_Graph  Movie . txt  , line_3 , ' - 
append ' ,  ' delimiter  1 ,  '  ' ) 

dlmwrite (' Line_Graph  Movie.txt  , movie  matrix,  - 

append ' , ' delimiter '  , ' \t ' , 'precision', ' % . 8f ' ) 
dlmwrite (' Line_Graph  Movie . txt  , line  end, 
append' , 'delimiter' , ' ' ) 

end 

end 

figure (2 ) 

plot (x  axis, S_plot, x^axis, S_f_plot, x^axis, S_conv_HT_plot, . . . 
x  axis , S_nideal  HT_plot,x  axis,S_ideal  HT_plot, . . . 
x  axis, S_HT_plot, x  axis, S_shock_plot, ' LineWidth ' , 3) 
title ('Total  Entropy  Flow  Rate  vs.  Axial 

Location ',' FontName ',' Times  New  Roman ',' Fontsize ', 24 ) 
legend (' Total  Entropy ',' Friction ',' Conv .  HT ',' Nonideal  HT ' , 'Ideal 
HT '  ,  ' Total  HT ' ,  ' Shocks ' ) 

xlabel (' Axial  Location  [m] ',' FontName ',' Times  New 
Roman ' , ' Fontsize ' , 20 ) 

ylabel (' Entropy  [ J/kg-K] ',' FontName ',' Times  New 
Roman ' , ' Fontsize ' , 20 ) 

set (gca, ' FontName ',' Times  New  Roman ',' Fontsize ', 20) 
figure  ( 3 ) 

plot (time, thrust_ss, time, thrust_t) 

legend (' Steady-State  Thrust ',' Transient  Thrust') 

title (' Thrust  vs.  Time ',' Fontsize ', 12 ) 

set (gca, ' FontName ',' Times  New  Roman ',' Fontsize ', 12 ) 

toe 

fprintf (' Creation  of  Line  Graph  Movie  Completed  \n  \n') 

end 


fprintf ( ' S  e  : 

%4 . 6f 

\n ' , S  plot (end) ) 

fprintf ( ' S  f  e  : 

%4 . 6f 

\n ! , S  f  plot (end) ) 

fprintf ( ' S  conv  HT  e: 

%4 . 6f 

\n',S  conv  HT  plot (end)) 

fprintf (' S  nideal  HT  e: 

%4 . 6f 

\n',S  nideal  HT  plot (end)) 

fprintf (' S  ideal  HT  e: 

%4 . 6f 

\n',S  ideal  HT  plot (end)) 

fprintf ( ' S  shock  e: 

%4 . 6f 

\n ' , S  shock  plot (end)) 

fprintf ( ' S  i : 

%4 . 6f 

\n ' , S_plot (1) ) 

fprintf ( ' S  before  bleed: 

%4 . 6f 

\n ' , S  plot  (400) ) 

fprintf (' S  after  bleed: 

%4 . 6f 

\n ' , S  plot  (501) ) 

fprintf ( ' S  lost: 

%4 . 6f 

\n',S  plot(400)-S  plot(501 

fprintf (' Percent  S  lost: 
S  plot(501))/S  plot(400) 

%4 . 6f 

) 

\n ' , (S  plot (400) - 

APPENDIX  C. 

ENTROPY  CALCULATIONS 
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Entropy  Generation  due  to  friction, 
1  r 


where, 


Sf=  j  [clx  xX(  (ut-u)+dx  xXb(u-ub)+dy  Tyr  (vr-v)+dy  xy|  (v-v,) 


du/ 


Entropy  Generation  due  to  convective  heat  transfer, 
1 


*con.  HT 


dy(k0T/ax)r-dJ'  (k  81  / 5x),+dx  (k  aT/ay)t  -dx  (k  aT/ey)b. 


where, 


k  klam+kturb. 


Entropy  Generation  due  to  ideal  heat  transfer, 


c  _  6Q 


Entropy  Generation  due  to  non-ideal  heat  transfer, 

SnHT=5Q  (™). 

Total  entropy, 

S  total  =  MFN2SN2  +  MFq2So2 

where, 

Sn2  =  SrefN2  +  aN2  ~  ^Ve/)  +  \CN2{j2  —  T'ref) 

+  jiN2(T3  -  T?ef)  +  y„2(r4  -  Tr%)  -  R„2  iog(£-) 
—  i?W2log  ( XNl ) 

^o2  =  Sreto2  +  ao2  log^r-^  +  bo2(T  —  Tre +  -Co2(T2  —  Tre 

+  \d02 ( T 3  -  T?ef)  +  \e02 (T*  -  T*„)  -  R0l  log 
~  «o2l°g  (*iv2) 
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Entropy  due  to  shocks, 


5shock 


S total  (■5/  "E  Sconv.HT  E 


Constants: 

aNl  —  3.1459RNz 
bN2  =  9.9154x10"4RN2 
cWz  =  -2.2912xlO~7RNz 
dN2  =  1.2181x10_11RN2 
eN2  =  1.1024x10-15RN2 
cl02  =  3.0809RNz 
b02  =  1.6962X10-3ROz 
Cq2  =  —  7.6334x10_7Ro2 
d02  =  1.7140xl(r10Ro2 
e0z  =  —  1.4116x10_14Ro2 
RNz  =  r/mwNz 
R0z  =  R/MWq2 

R  =  8,314.34  Vkmol  K 

MWNz  =  28.016  kgNz/kmolN2 

MW0z  =  32.000  kgo2/kmol02 

MFNz  =  0.7686  kgNz/kgair 

MF02  =  0.2314  kg0z/kgair 

XNz  =  kmolesN2/(kmolesNz  +  kmoles02) 

Xq2  =  kmoleso2/(kmoleso2  E  kmoles02) 

SrefN2  =  1.9161x105/MWNz 

Sref02  =  2.05142x105/MWOz 

Pref  =  100,000  pascals 

Tref  =  298.0  Kelvin 


APPENDIX  D. 

QUASI-ONE-DIMENSION  ANALYTICAL  THROTTLING  MODEL 
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function  min  max_throttle 

clear  all 
close  all 
clc 

%Given 

gamma  =  1.4; 

R  =  287; 

C_p  =  gamma/ (gamma  -  1)  *R; 

M_0  =  3.0; 

P_0  =  26436.267594; 

T_0  =  223.15; 

rho_0  =  P_0/(R*T_0); 

u  0  =  M  0*sqrt (gamma*R*T  0); 

h~=  4 . 4  2  e  7 ; 

r_i  =  0.0625/4; 

A  =  pi*r_iA2; 

AR_1  =  0.6227; 

AR_2  =  0.7464; 

A_m  =  AR_1 *A; 

A_t  =  AR_2*A; 
mdot  =  rho  0*u  0*A 

T_t_0  =  T_0* (1  +  (gamma  -  1)/2*M_0A2); 

P  t__0  =  P  0*(1  +  (gamma  -  1)/2*M  0A2 ) A (gamma/ (gamma  -  1)); 

%--  Start  calculations  for  max  thrust  (normal  shock  at  A  m)  - 

%First  inlet 

A  A  star  bs  =  ( (gamma  +  1) / 2)  A (- (gamma  +  1) / (2* (gamma  -  1) ) ) * (1  + 

(gamma  -  1) /2*M^0A2) A ( (gamma  +  1) / (2* (gamma  -  1) ) ) /M_0; 

A  star  bs  =  A*l/A  A  star  bs; 

A  m  A  star  bs  =  A  m/A  star  bs; 

f  M  m  bs  =  @ (M  m  bs)  ( (gamma  +  1) /2) A (- (gamma  +  1) / (2* (gamma  -  1) ) ) * (1 
+  (gamma  -  1)/2*M  m  bsA2)A( (gamma  +  1) / (2* (gamma  - 
1) ) ) /M  m  bs  -Am  A_star_bs; 

M  m  bs  =  fzero(f  M  m  bs,2.2610); 

P  t  as  =  P  t_0* (  ( (gamma  +  1) / 2*M  m  bsA2) / (1  +  (gamma  -  1) /2*M  m  bsA2) 

) A (gamma/ (gamma  -  1) ) * (  1/ (2*gamma/ (gamma  +  1) *M  m  bsA2  - 
(gamma  -  1)/ (gamma  +  1))  ) A (1/ (gamma  -  1)); 

M  m  as  =  sqrt (  (M  m  bsA2* (gamma  -  1)  +2)/ (2*gamma*M  m  bsA2  -  (gamma  - 

l ) )  )  ; 

A  m  A  star  as  =  (  (gamma  +  1)  / 2)  A  (-  (gamma  +  1)  /  (2*  (gamma  -  1)  )  )  *  (1  + 
(gamma  -  1)/2*M  m  asA2)A( (gamma  +  1) / (2* (gamma  - 
1 ) ) ) /M_m_as; 

A  star^as  =  A  m*l/A  m  A  star  as; 

A_A_star_as  =  A/A_star_as; 

f  M  2  =  @  (M  2)  (  (gamma  +  1)  / 2)  A  (-  (gamma  +  1)  /  (2*  (gamma  -  1)  )  )  *  (1  + 

(gamma  -  1) /2*M  2A2) A ( (gamma  +  1) / (2* (gamma  -  1) ) ) /M  2  - 
A_A_s  t  a  r_a  s ; 

M  2  =  fzero(f  M  2,0.2352); 

M_3  =  M_2; 

P_t_3  =  P_t_as; 

P  3  =  P  t  3/ (1  +  (gamma  -  1) /2*M  3A2) A (gamma/ (gamma  -  1) ) ; 
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%Second  inlet  and  Rayleigh  burner 


f  M  4  =  @ (M  4)  ( (gamma  +  1) /2) A (- (gamma  +  1) / (2* (gamma  -  1) ) ) * (1  + 

(gamma  -  1) / 2*M  4A2) A ( (gamma  +  1) / (2* (gamma  -  1) ) ) /M  4  -  A/A  t; 
M  4  =  fzero(f  M  4,0.5034); 

P  4  =  P_3*(l  +  gamma*M_3 A2 ) / ( 1  +  gamma*M  4A2); 

Pt4=P4*(l+  (gamma  -  1)/2*M  4A2)A (gamma/ (gamma  -  1) ) ; 
f_2  =  (  (1  +  (gamma  -  1) /2*M_2A2) / (1  +  gamma*M_2 A2 ) A2  )*M_2A2; 

f_4  =  (  (1  +  (gamma  -  1) /2*M_4A2) / (1  +  gamma*M_4 A2 ) A2  )*M_4A2; 

T_t_4  =  T_t_0  *  f _4 / f _2 ; 

f_M_e  =  @ (M_e)  ((gamma  +  1 ) /2 ) A (- (gamma  +  1)/ (2* (gamma  -  1)))*(1  + 

(gamma  -  1) /2*M_eA2) A ( (gamma  +  1) / (2* (gamma  -  1)))/M  e  -  A/ A  t; 
M_e  =  f zero (f_M_e, 1 . 6955) ; 

P  e  =  P  t  4/ (  (1+  (gamma  -  1) /2*M  eA2) A (gamma/ (gamma  -  1) )  ) ; 

T_e  =  T_t_4/ (  (1  +  (gamma  -  l)/2*M_eA2)  ); 

u  e  =  M  e*sqrt (gamma*R*T  e)  ; 

F  max  =  mdot* (u  e  -  u_0 )  +  (P  e  -  P_0)*A; 
mdot_f_max  =  mdot/h* ( f_4 / f _2  -  l)*C_p*T  t  0; 

Qdot  max  =  mdot  f  max*h; 


% -  Start  calculations  for  min  thrust  (normal  shock  at  A_2 )  - 

%First  inlet 
M_2Jos  =  M_0 ; 

P  t  as  =  P  t_0* (  ((gamma  +  1)/2*M  2_bsA2)/(l  +  (gamma  -  1)/2*M_2  bsA2) 

) A (gamma/ (gamma  -  1) ) * (  1/ (2*gamma/ (gamma  +  1) *M  2  bsA2  - 
(gamma  -  1)/ (gamma  +  1))  ) A (1/ (gamma  -  1)); 

M  2  as  =  sqrt (  (M  2  bsA2* (gamma  -  1)  +  2 ) / (2*gamma*M  2  bsA2  -  (gamma  - 
1)  ) 

M_3  =  M_2_as ; 

P_t_3  =  P_t_as; 

P3=Pt3/(l+  (gamma  -  1) /2*M  3A2)A (gamma/ (gamma  -  1) ) ; 


%Second  inlet  and  Rayleigh  burner 

f  M  4  =  @ (M  4)  ( (gamma  +  1) / 2)  A (- (gamma  +  1) / (2* (gamma  -  1) ) ) * (1  + 

(gamma  -  1) /2*M  4A2) A ( (gamma  +  1) / (2* (gamma  -  1) ) ) /M  4  -  A/A  t; 
M  4  =  fzero(f  M  4,0.5034); 

P  4  =  P_3*(l  +  gamma*M_3 A2 ) / ( 1  +  gamma*M  4A2); 

Pt4=P4*(l+  (gamma  -  1) /2*M  4A2)A (gamma/ (gamma  -  1) ) ; 
f  2  =  (  (1  +  (gamma  -  1)/2*M  2_asA2)/(l  +  gamma*M_2  asA2)A2  ) *M_2  asA2; 
f_4  =  (  (1  +  (gamma  -  1) /2*M^4A2) / (1  +  gamma*M_4 A2 ) A2  )*M_4A2; 

T_t_4  =  T_t_0  *  f _4 / f _2 ; 

f_M_e  =  @ (M_e)  ((gamma  +  1 ) /2 ) A (- (gamma  +  1)/ (2* (gamma  -  1)))*(1  + 

(gamma  -  1)/2*M  eA2)A(  (gamma  +  1)/ (2*  (gamma  -  1)))/M  e  -  A/A_^t; 
M_e  =  f zero (f_M_e, 1 . 6955) ; 

P  e  =  P  t  4/ (  (1+  (gamma  -  1) /2*M  eA2) A (gamma/ (gamma  -  1) )  ) ; 

T_e  =  T_t_4/ (  (1  +  (gamma  -  l)/2*M_eA2)  ); 
u  e  =  M  e*sqrt (gamma*R*T  e)  ; 

F  min  =  mdot* (u  e  -  u_0 )  +  (P  e  -  P_0)*A; 
mdot_f_min  =  mdot/h* ( f_4 / f _2  -  l)*C_p*T  t  0; 

Qdot  min  =  mdot  f  min*h; 
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